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ABSTRACT 


AGGIE I is a computer program for predicting the linear and nonlinear, 
static and dynamic structural response of two- and three-dimensional 
continuum solids. The program is based on, isoparametric finite elements 
and allows for 2-D plane stress, plane strain, and axisymmetric analyses 
and general 3-D analyses. Large strain kinematics is based on the total 
Lagrangran formulation. Materially nonlinear models include several elastic- 
plastic work-hardening models as well as an incompressible Mooney-Rivlin model. 

Included in this report is a brief description of the theoretical bases 
of the program, the material models used, the element library and the overall 
program organization. Instructions for data input preparation are given in 
detail. Several sample problems are given along with the required program 
input and program generated solutions. 
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1. INTRODUCTION 


The development of computer programs for nonlinear structural analyses 
has progressed at a steady pace for the past ten years beginning with 
special application programs and progressing to larger, more general 
purpose programs like NONSAP, ANSYS, NASTRAN, and MARC. 1 Many of the 
larger general purpose programs which are currently available were originally 
developed for linear structures and later extended to include a nonlinear 
capability. The research reported herein was aimed at developing a mod- 
erately sized finite element program specifically for nonlinear structural 
analysis. 

In beginning the development of such a program, a review of existing 

nonlinear finite element programs was undertaken with the purpose being 

to determine what analysis capability currently existed, what techniques 

were being used, and to determine if any existing computer program might 

be suitable as a beginning base program. From this review, it was felt 
2 

that the NONSAP program developed by Bathe, Wilson and co-workers would 
form a suitable base from which to begin the work. The NONSAP program 
was modularized making it easy to modify, the program organization and 
logic was good, it was fairly well documented and it had a good library 
of material models and elements to conduct geometrically (large strain) 
and materially nonlinear analyses of two-dimensional structures (the 3-D 
capability was limited to linear elastic problems only). 

Using the NONSAP program as a base, a moderately sized finite element 
program for nonlinear structural analysis has been developed. The pro- 
gram is based on two- and three-dimensional isoparametric solid elements. 
Nonlinearities accounted for in the program may be due to large dis- 
placements, large strains and nonlinear material behavior (linear be- 
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havior is also allowed). A number of material models have been in- 
corporated; these include linear elastic, orthotropic, nonlinear elastic, 
nonlinear elastic incompressible and elastic-plastic models. Both 
static and dynamic response predictions are possible. The program also 
allows for the possibility of displacement incrementation for use in 
predicting buckling and post-buckling response. Dynamic dimensioning 
of most FORTRAN arrays allows efficient utilization of available high- 
speed memory. The program uses either an out-of-core or an in-core 
equation solver depending upon problem size and available main memory. 

The mesh description (node and element information) is input on a node 
by node and element by element basis or may be generated by using a 
rudimentary substructure option. Provision is made for both conservative 
and nonconservative (deformation dependent) loading. 

A brief description of the theoretical basis of the program, the 
elements used, the available material models, and program organization is 
given in the next section. More detailed information regarding these 
items may be found in Refs. 2-5. A description of the input data required 
by the AGGIE I program is then presented. Sections are then presented 
which give suggestions and guidelines for use of the program#- program 
resource requirements, capacity, limitations, and the program' restart 
option. A last section contains selected example problems along. with 
listings of the required input data and program results. 
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2. SUMMARY OF AGGIE I COMPUTER PROGRAM 

The purpose of this section is to provide a summary of the theoret- 
ical basis of the computer program, the elements used, the available 
material models, and the program organization. Details of these may be 
found in Refs. 2-5. 


Incremental Equilibrium Equations 

The linearized incremental equilibrium equations which are solved 
at each time or load step are given by 

[M]{u t+At } + [C]{ii t+At } + P^Hau} = {R t+At } - {F t } (1) 

where 


[M] = mass matrix 
[C] = damping matrix 

4 - 

[K ] = tangent stiffness matrix at time t 

{pt+Atj _ externa ] noc | a ] f 0rce vector at time t + At 

t 

{F } = nodal force vector due to internal element stresses at time t 
{u > = nodal acceleration vector at time t + At 
{au} = nodal displacement increment vector defined by u t+At = u*" + au 
For a static problem, the mass and damping matrices are omitted and the 
time increment At may be interpretted as an equivalent load increment. 

The incremental equations given above yield only an approximate 
solution (because of the linearization used in their derivation) and 
hence will yield solutions which have accumulated error after several 
steps. To ensure that equilibrium is being satisfied at each step, 
one may use a modified Newton-Raphson method to perform equilibrium 
iteration. The iterative equations take the form 

[M]{ii t+At }^ ^ + [C]{u t+At } (l) + [K^Iaau}^ = {R t+ At}(i ~1 ) - {F t+At } (i-l) ^ 
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where (i) denotes the i th equilibrium iteration and {aau}^ denotes 
the correction to the displacement increment at iteration i, that is, 

{Au}( 1+1 ^= {Au } 1 + {aau } 1 


{ u t+At }i +1 


= {U t } 1+1 + {AU } 1+1 


It should be pointed out that in a general Newton type iteration, the 
tangent stiffness matrix [K ] would be updated at each iteration 
(i.e., use [K t+At ]^); however, in the modified Newton method, this 
term is held constant during the iteration (or updated only periodically 
when convergence deteriorates). Input parameters allow the program user 
to specify the updating frequency. 

The time integration procedure used in dynamic analyses is either the 
Wilson e or the Newmark g methods. 

The original NONSAP program utilized both the total Lagrangian anct 
updated Lagrangian formulations in developing the incremental equilibrium 
equations. All element and material models developed during the current 
research program have been restricted to the total Lagrangian formulation; 
however, those material models in the original NONSAP program which 
permitted the use of the updated Lagrangian solution have been retained 
and are still available in AGGIE I. 

Element Library 

Several types of elements are currently available in the AGGIE I 
program. These include a three-dimensional truss element, two-dimensional 
isoparametric solid element with 4 to 8 nodes, and a three-dimensional 
isoparametric solid element with 8 to 21 nodes. The two-dimensional 
element is reducable to a constant strain triangular element and the 
three-dimensional element is reducable to a constant strain tetrahedron. 
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The tfuss element is shown in Fig. 1 and has in general three degrees- 
of- freedom per node. Material properties and cross-sectional area are 
assumed to be constant over the length of any element. The element may 
be used in linear elastic or nonlinear elastic analyses in combination 
with small or large displacements (for a large displacement analysis, only 
the updated Lagrangian formulation is available). 

The general two-dimensional isoparametric solid element is shown in 
Fig. 2. The element may be used for plane stress, plane strain or axi- 
symmetric analyses. The material models which are available for the two- 
dimensional element are shown in Fig. 2 and are described briefly in the 
next section. This element may be used to model axisymmetric shells; 
however, because of the continuum formulation used, the element does not 
represent Kirchh off- type, thin shell behavior very well. 

The general three-dimensional isoparametric solid element 
is shown in Fig. 3. The element may have from 8 to 21 nodes with 3 
degrees -of-freedom per node. The material models which are^currently 
available for the three-dimensional element are also listed in Fig. 3. 

Material Models 

The original NONSAP program contained a number of material models 
for the two-dimensional elements but contained only a linear and nonlinear 
elastic material model for the three-dimensional elements. In the AGGIE 
program, new material models have been added for the two-dimensional 
elements and four material models have been added for the three-dimensional 
elements. Table 1'gives a complete listing of the available material models. 

For the truss elements , two- models are available: 1) linear elastic 

material and 2) nonlinear elastic material. The nonlinear elastic-material 
model makes use of a user-specified, piece-wise linear stress-strain curve. 
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Geometric Nonlinear 

Formulations 

Available Material Models 

Li near 

— 

■SRSjwBSI 

Updated 

Lanqranqian 

1. linear elastic 

2. nonlinear elastic 

X X 

1 


X 

X 


Fig. 1 Three-Dimensional Truss Element 
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Geometric Nonlinear Formulation 


Available Material Models 


Total 

1 

Updated 

Linear 

Laqranqian 

Laqranqian 

1 . 

linear elastic, isotropic 

X 

X 

X 

2. 

linear elastic, orthotropic 

X 

X 


3. 

variable tangent moduli model 

X 

X 


4. 

curve desicription model (plane 
strain, axi symmetric) 

X 



5. 

curve description model with 
tension cutoff 

X 



6. 

elastic plastic, combined 
ki nemati c-i sotropi c hardeni ng , 
von-Mises yield condition 

X 

X 


7. 

elastic-plastic, Drucker-Prager 
yield condition 

X 



8. 

incompressible, nonlinear elastic, 
Mooney-Rivlin material (plane 
stress) 

X 

X 


9. 

incompressible, nonlinear elastic, 
Mooney-Rivlin material (plane 
strain or axi symmetric) 

X 

X 


10. 

elastic-plastic, mechanical 
sublayer model, von-Mises yield 
condition 

X 

X 



Fig. 2 Two-Dimensional Isoparametric Continuum Element for 
Plane Stress, Plane Strain or Axi symmetri c Analysis 
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Geometric Nonlinear Formulation 

Available Material Models 

Linear 

Total 

Laqranqian 

Updated 

Laqranqian 

1. linear elastic, isotropic 

X 

X 


2. linear elastic, orthotropic 

X 

X 


3. curve description model 

X 

X 


4. elastic-plastic, combined 

X 

X 


ki nemati c-isotropi c hardening 
von-Mises yield condition 




5. elastic-plastic, mechanical 

X 

X 


sublayer model, von-Mises 
yield condition 




6. incompressible nonlinear 

elastic, Mooney-Ri vlin material 

X 

X 



Fig. 3 Three-Dimensional Isoparametric Continuum Element 
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Table 1. Material Models 


Model 8 1 -D 


1 linear 

elastic 


2 | nonlinear 

f elastic 



2-D 

3-D 

linear elastic, 
isotropic 

linear elastic, 
isotropic 

linear elastic, 
orthotropic 

linear elastic, 
orthotropic 

variable tangent 
moduli model 

curve description 
model 

curve description 
model (plane strain, 
axi symmetric) 

.elastic-plastic, 
combined kinematic- 
isotropic hardening, 
VonMises yield 
condition 

curve description 
model with ten- 
sion cutoff 

elastic-plastic, 
mechanical sub- 
layer model , Von- 
Mises yield condition 

elastic-plastic, com- 
bined kinematic- 
isotropic hardening, 
Von Mises yield 
condition 

incompressible nonlinea! 
elastic, Mooney-Rivlii 
material 

elastic-plastic, Drucker 

Prager yield conditio 



n 

i 

incompressible, nonlinear 
elastic, Mooney-Rivlin 
material (plane stress) 

incompressible, nonlinea 
'elastic Mooney-Rivlin 
material (plane strai 
! or axi symmetric) 

r 



1 

elastic-plastic, mechanical 
sublayer model, Von Mises 
yield condition | 
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For the two-dimensional solid elements , a number of models are available 
as shown in Table 1. Most of the material models may be used with either a 
geometrically linear (small strain) or nonlinear (large strain) formulation. 

Two linear material models are available, i.e., either isotropic or 
orthotropic. 

A newly developed nonlinear elastic (hyperelastic) incompressible 

model which makes use of Mooney-Rivlin material constants to describe the 

nonlinear stress-strain relation may be used for the analysis of some 

rubber-like materials. The plane-strain and axisymmetric models make use 

of Lagrange multiplier constraints in order to impose the incompressibility 
5 

condition. 

Two new elastic-plastic material models provide a wide range of 
hardening rules which may be used in elastic-plastic analyses. Both 
models are based on the von Mises yield condition and the associated 
flow rule. For the hardening rule, one model makes use of a combined 
kinematic-isotropic hardening model while the second is based on a me- 
chanical sublayer (Besseling) model. The combined kinematic-isotropic 
hardening model allows for the input of a yield surface size vs. equiv- 
alent uniaxial plastic strain curve if such data is available; otherwise 
this information is generated by the program for a user-specified value 
of the ratio of kinematic-to-isotropic hardening. Both models make use 
of the von Mises yield condition and associated flow rule. Details of 
both models may be found in Ref. 4. 

Three geological material models originally in NONSAP have been re- 
tained in the present program. These are the variable tangent moduli 

model, the curve description model and the curve description model with 
3 

tension cutoff. The variable tangent moduli model describes an isotropic 
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material in which the bulk and shear moduli are functions of the stress 
and strain invariants and the functional relationship used replaces an 
explicit yield conditions. In the curve description model, the instant- 
aneous bulk and shear moduli are defined by piecewise linear functions 
of the current volume strain. With the tension cutoff model, tensile 
stresses due to applied loading cannot exceed the gravity in-situ pressure; 
the model assumes reduced stiffenss in the direction of the tensile stress 
which exceed in magnitude the gravity pressure (i.e., when a "crack" 
develops) . 

The material models available for the three-dimensional solid elements 
are also shown in Table 1. Several new models were developed and implemented 
in the base program. These include the Mooney-Rivlin incompressible model, 
the linear orthotropic model, the elastic-plastic model with combined 
kinematic-isotropic hardening, and the elastic-plastic model with mechanical 
sublayer models. Each of these models may be used in either a geometrically 
linear analysis or a large strain analysis based on the total Lagrangian 
formulati on. 

Program Organization 

The general program organization which was utilized in the base program 
NONSAP^ has been retained in the present program. The significant dif- 
ferences are in the way the element matrices are assembled into the 
structural matrices and the allocation of the high-speed storage. 

As in NONSAP, the solution process is divided into three phases: 

!• Input and Mesh Generation Phase . Various control information 
is first read from data. Depending upon which input option is selected, 
the node and element information is either read from data or generated 
internally and then stored on tape or disk. Additionally, equation numbers 



for all active degrees of freedom are generated and element-node con- 
nectivity arrays are determined and stored on disk. Externally applied 
pressure loads are converted to nodal forces and stored on disk. 

2. Assemblage of Constant Structural Matrices . The stiffness, 
mass and damping matrices for all linear elements are computed, assembled 
and stored on disk. 

3. Step-by-Step Solution . The incremental equations of equilibrium 
are solved for each step. At each step, equilibrium iteration is per- 
formed (if desired) until a converged solution is obtained. Displacements, 
velocities, accelerations, stresses and strains are computed and printed. 

The program is designed to process elements by groups; all elements 
within an element group have like characteristics, i.e., they are all the 
same element type (e.g., 2-D axisymmetric) , they all have the same type 
of geometric nonlinearity and the same material model. Each element group 
generally calls for a specific overlay segment of the program. Program 
efficiency is increased by grouping elements to avoid unnecessary or 
repeated calling of overlay segments. Element groups may be either linear 
(geometrically and materially linear) or nonlinear. During all phases, 
linear element groups are processed first. Stiffness properties for all 
linear element groups are computed during phase 2 and stored on disk. 

During phase 3, the structural stiffness for all linear element groups 
is read from disk; then stiffness matrices for all nonlinear element groups 
are computed (based on the current stress strain-state) and added to the 
linear contribution. 

The program is designed with the idea of making maximum use of avail- 
able high speed memory. A common block /AMAT/ dimensioned for MTOT lo- 
cations is used to store the arrays that are required during the different 
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phases of program execution. During the input phase, space (NUMEST locations) 
is allocated in the common block for storage of element property data 
(element coordinates, element-node connectivity, material properties and 
working arrays for element stresses, strains and various variables asso- 
ciated with nonlinear material models). For each element group, the input 
is processed using the same NUMEST locations and necessary information 
is written on disk. During the solution phase, the element property data 
is read from disk and the same space is used repeatedly as each element 
group is processed. The program automatically allocates the NUMEST lo- 
cations within the common block based on the storage requirements of the 
largest element group. During the input phase, the program determines 
the storage that will be required for the structural coefficient (e.g., 
stiffness in static analyses) matrix and the space (LBLOCK) that will be 
available with the common block. If sufficient space is available in the 
common block to store the coefficient matrix in core, then this is done. 

If sufficient space is not available, then the coefficient matrix is 
stored on disk in blocks of length LBLOCK and then processed block-by- 
block within the LBLOCK locations of common block AMAT. More detailed 
information on storage allocation may be found in the Programmer's Guide 
for the program. 

As mentioned before, the AGGIE program contains over 125 subroutines 
and is overlaid. The overlay structure currently v being used is shown in 
Table 2 (although additional overlaying may be possible and desirable). 

Equation Solvers 

The solution of the incremental equilibrium equations is obtained 
using either the in-core linear equation solver 0PTS0L 7 or the out-of- 

O 

core linear equation solver OPTBLK. The coefficient matrix (left side 
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Table 2 

Overlay Structure of AGGIE I Program 













































of the equations) is stored as a one-dimensional array using a skyline 
storage scheme (by columns) i.e., only the terms in each column which 
fall below the skyline are stored. If the coefficient matrix will fit 
within the high-speed memory space which has been allocated to it, then 
the solution of the equations is done in-core by OPTSOL. If the coef- 
ficient matrix will not fit entirely in-core, then it is stored on disk, 
tape or other low-speed storage device by blocks and solved by QPTBLK 
in a block-by-block fashion. The blocksize is selected as large as pos- 
sible depending upon available high-speed memory. 
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3. FUTURE IMPROVEMENTS 


The AGGIE program is currently being extended to include the 
following improvements: 

1. General 3-D thin shell element 

2. General 3-D beam element 

3. Thermal strains 

4. Non-isothermal , combined kinematic-isotropic hardening 
model with temperature dependent material properties 

5. Creep models 

6. Nonlinear viscoelasticity 

Additional improvements to be added in the future include: improved 

mesh generation, non-zero boundary conditions, and additional element 
geometries . 

Users interested in obtaining further information on the status 
of current program improvements should contact: Dr. Walter E. Haisler, 

Aerospace Engineering Department, Texas A&M University, College Station, 
Texas 77843, Phone (713) 845-7541. 
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4. PROGRAM INPUT 


This section describes the required input data for the AGGIE program. 
The input data is separated into 16 sections; some of these are generally 
not required for a particular problem. Multiple cases may be run by 
inputting the required input sequentially for each case. The data set 
for the last problem case should be followed by two blank cards; however, 
omission of these will only cause an end-of-file read error when the 
program tries to process the non-existent data set. An error in one 
input data set will cause the program to terminate and not process any 
data sets which might follow. 

The input instructions here should be used in conjunction with the 
example problems and deck setups shown in a later section. Additionally, 
a section follows which outlines some of the known limitations of the 
program and gives some suggestions which may be useful when difficulties 
are encountered in using the program. 

The input data consists of the following sections: 

I. Heading Card 

II. Program Control Cards 

III. Nodal Point Data 

IV. Truss Elements 

V. Two-Dimensional Solid Elements 

VI. Three-Dimensional Solid Elements 

VII. Node and Element Input by Substructures 

VIII. Load Curves 

IX. Nondonservati ve Pressure Loads 

X. Conservative Pressure Loads 

XI. Concentrated Nodal Forces 

XII. Rayleigh Damping 

XIII. Concentrated Nodal Masses 

XIV. Concentrated Nodal Dampers 

XV. Initial Conditions 

XVI. Frequency Solution Data Card 
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As noted previously, the AGGIE program is an extension of the NONSAP 
program. Hence, some sections of the input described here are necessarily 
the same as described in the NONSAP user's manual (with only minor mod- 
ifications). This is the case for input instructions which refer to 
program segments and material models which were originally in NONSAP and 
were left intact and whose input instructions were unchanged*, e.g. , 
nodal point input, all linear and geological material models, initial 
conditions, and concentrated mass and damping input. These input sections 
are marked with an asterisk (*). 


18 



I. HEADING CARD (20A4) 


column variable description notes 

1-80 HEAD (20) Heading information which will (1) 

appear on output. 

NOTES 

(1) Any information placed on this card will appear on the title page 
of the output. 
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II. PROGRAM CONTROL CARDS 

II. 1 Solution Control Variables (611, 14, 415, 2F1Q.0, 315, 612) 


columns 

variable 

description 

notes 

1 

ID0F(1) 

Master X-translation code 
EQ.O - admissable 
EQ.l - deleted 

(1) 

2 

ID0F(2) 

Master Y-translation code 


3 

ID0F(3) 

Master Z-translation code 


4 

ID0F(4) 

Master X-rotation code 


5 

I D0F( 5) 

Master Y-rotation code 


6 

I D0F( 6) 

Master Z- rotation code 


7-10 

NUMEG 

Number of substructures or element groups 

(2) 

11-15 

MGO 

Mesh generation option 
EQ.O - Individual node and element data 
input 

EQ.l - Mesh described by substructures 

(3) 

16-20 

MODEX 

Solution mode 
EQ.O - Data check only 
EQ.l - Execution 
EQ.2 - Restart 

(4) 

21-25 

IDINC 

Displacement incrementation Flag 

EQ.O - Load incrementation 

EQ.l - Displacement incrementation 

(5) 

26-30 

NSTE 

Number of solution steps 

(6) 

31-40 

DT 

Load (time) step increment 

'6) 

41-50 

TSTART 

Load (time) value at which solution 
starts 

(7) 

51-55 

IMASS 

Flag indicating static or dynamic analysis (8) 

EQ.O - static analysis 

EQ.l - dynamic analysis with lumped, diagonal mass 
EQ.2 - dynamic analysis with consistent mass matrix 

56-60 

IEIG 

Frequency solution flag 
EQ.O - no frequency solution 
EQ.l - perform frequency analysis 

(9) 

61-65 

IPLOT 

Plot tape flag 

EQ.O - no plot tape 

EQ.l - plot tape generated 

(10) 
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II. PROGRAM CONTROL CARDS (continued) 


columns variable 


description 


notes 


66-67 

68- 69 

69- 70 

70- 71 

71- 72 
76-77 


ICOMB(l) Substructure combining flags. Used with 
ICOMB (2) substructure mesh generation only (MG0=1). 
ICOMB ( 3) 

ICOMB (4) 

ICOMB (5) 

ICOMB (6) EQ.O - Program generates one element group 
for every substructure. 

GT.O - Program combines substructures 
IC0MB[I(odd)] to ICOMB [1+1 (even)] into one 
element group. Three sets of combinations 
are allowed. Leave blank if no substruc- 
ture combining is desired. 


NOTES 

(1) The master translation codes are used to delete nodal degrees-of- 
freedom (dof) which are not required for this problem. The 
program allows for six degrees of freedom per node; however, since 
many problems do not require the complete set, these master 
translation codes are used to delete these dof and thereby con- 
serving storage and- execution time. Any later referenc to a 
deleted degree-of-freedom is ignored by the program. 

All elements currently available in this program allow for 
translational degrees-of- freedom only and; hence, IDOF (I) 
for 1=4,5 and 6 must equal 1. The 2-D continuum element is 
assumed to lie in the (Y-Z) plane and; hence, I DOF (1 ) must equal 1 
for this case. 

(2) The input of element information is handled by element groups 
each of which has like characteristics as discussed later. If the 
nodal mesh is described by substructures (MG0=1), then NUMEG 
represents the number of substructures by which the complete 
structure geometry is described and which are read in Section VII. 

If the nodal mesh is input by individual nodes, then NUMEG represents 
the number of element groups to be input in Sections IV - VI. 

In either case, the elements comprising an element group or sub- 
structure must have the same characteristics , i.e., element type, 
type of analysis, material model and numerical integration order 
(these are described by the NPAR card in Section V.l or VI. 1). 

(3) The nodal coordinates and element-node connectivity may be 
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II. PROGRAM CONTROL CARDS (continued) 


described by two approaches. The node numbers and coordinates 
may be input individually as described in Section III. Al- 
ternately, the geometry and element-node connectivity may be 
generated by describing “substructures" and letting the program 
automatically generate internal node numbers, coordinates, element 
numbers and element-node connectivity information. This approach 
is described in Section VII. 

(4) In the data check mode, all input is read, partially checked for 
validity, and printed. The restart mode is used to restart a 
problem for which a prior solution was obtained and a restart 
tape written. Detailed instruction on restart jobs is contained 
in Section PROGRAM RESTART. 

(5) For the ususal case, load (time) is considered to be the inde- 
pendent variable and, hence, the load is incremented from zero 
to the final value. However, in static buckling problems it 
may be more convenient to let some characteristics displacement 
be the independent variable and increment this displacement 
(and solve for the corresponding load). 

(6) For dynamic problems, the time step DT must be selected to 
maintain accuracy and stability. The total time for which a 
solution is obtained is given by NSTE*DT. The loads must be 
specified over this time range and are described by force vs. 
time tables input in Section VIII. 

For static problems the load increment DT is used to increment 
a load parameter. The loads are similarly described in terms 
of force vs. load parameter tables. 

For load incrementation, the load (time) increment DT is not used 
since the incrementation is controlled by the displacement 
increment DINC input in Section II. 7. 

(7) Normally, TSTART is set to zero. However, when restarting a 
problem, TSTART is set to the value of the load parameter (or 
time) at which the solution is to restart and for which restart 
information is available on tape from a previous run. 

(8) For dynamic analyses, inertia effects are accounted for and, in 
addition, damping may be included. 

(9) Mode shapes and frequencies are computed prior to the dynamic 
solution if I El G = 1 (Note - this part of the program is not 
operational in this version of AGGIE). 

(10) The program prepares a plot tape of initial mesh, geometry. Subroutine 
MESH2D utilizes the plot tape in combination vnth several TAMII Cal comp 
library programs to plot two-dimensional meshes. Although, a plQt tape 
is created for three-dimensional meshes, no plotting routines are 
included. Users utilizing the plot tape will undoubtedly have to 
interface the program with their installations plotting and graphics 
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II. PROGRAM CONTROL CARDS (continued) 


II. 2 Nonlinear Solution Variables (415, F10.0)* 


col umns 

variable 

description 

notes 

1-5 

ISREF 

Number of load steps between re- 
formation of stiffness matrix. 

If EQ.O - default reset to 1. 

0) 

6-10 

NUMREF 

Number of allowable stiffness re- 
formations in each step. 

EQ.O ■= default reset to 1 

(2) 

11-15 

IEQUIT 

Number of steps between equilibrium 
i terati ons . 

EQ.O - default reset to 1 

(3) 

16-20 

ITEMAX 

Maximum number of equilibrium 
iterations permitted before re- 
forming stiffness 
EQ.O - default reset to 15. 

(2) 

21-30 

RTOL 

Tolerance used to measure 
equilibrium convergence. 

EQ.O - default reset to 0.001 

(4) 


NOTES 

(1) For nonlinear element groups (materially or geometrically nonlinear), 
the stiffness matrix is reformed every ISREF load (time) steps based 
on information at the end of the previous step. In addition, 
history dependent variables are updated for use in any history 
dependent material law. 

(2) After ITEMAX equilibrium iterations in any load step, the stiffness 
matrix is reformed if NUMREF. GE.2 and iteration continued. The 
stiffness is always reformed at the beginning of the step prior to 
the first iteration (for those steps where reformation is to be 
performed) and then reformed NUMREF times with ITEMAX iterations 
between each reformation. 

(3) For nonlinear problems, equilibrium iteration is performed every 
IEQUIT steps until convergence of the displacements is obtained. 

For example, if IEQUIT is set to 3, then iteration is performed 
during load steps 3, 6, 9, etc. only until the solution converges to 
a tolerance of RTOL. If convergence does not occur within ITEMAX 
iterations, then the stiffness is reformed (if NUMREF. GT.l) and 
iteration continued; however, if stiffness reformation is not allowed 
then the solution is stopped. 


*An asterisk is used to mark those sections of the input which are 
identical to the NONSAP Program 



II. PROGRAM CONTROL CARDS (continued) 

(4) During equilibrium iteration, convergence of the solution is 

determined by the change in the Euclidean norm of the displacement 
vector, i .e. , when 


,(n-l) 






<RT0L 


where || u^ || is the Euclidean norm of the displacement vector 
at iteration "n". 
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II. PROGRAM CONTROL CARDS (continued) 


II. 3 Load Control Variables (515) 


col umns 

variable 

description 

notes 

1-5 

NLCUR 

Number of load curves to be input 
EQ.O - Linear load curve in internally 
generated from load value TSTART 
to TSTART + DT*NSTE . 

GT.O - Number of load curves to be 
read from data 

(1) 

6-10 

NPRESC 

Number of conservative pressure load 
cards to be input 

(2) 

11-15 

NPRESN 

Number of nonconservative pressure 
load cards to be input 

(2) 

16-20 

NCONCF 

Number of concentrated nodal force 
cards to be input 

(2) 

21-25 

NBODY 

Body force flag 

EQ.O - no body forces 

EQ.l - body forces to be included 

(2-D only) 

(3) 


NOTES 

(1) The load curves describe the variation of the load functions with 
load parameter (time) and are input in Section VII I . 

For displacement incrementation (IDINC.EQ.l ) , no load curves are 
input (set NLCUR.EQ.O). In displacement incrementation, the load 
curve multipliers are computed as part of the solution. 

(2) Conservative and nonconservative (deformation dependent) pressure 
loads applied to the face of an element or substructure and con- 
centrated nodal forces are input in Sections IX - XI. 

Nonconservative loads are not allowed for displacement incrementation 
(IDNIC.EQ.l). 

(3) Body forces may be applied to 2-D elements and are assumed to act 
in the global 2 direction. 
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II. PROGRAM CONTROL CARDS (continued) 


II. 4 Print Control Variables (815) 


col umns 

variable 

description 

notes 

1-5 

IPRI 

Output printing interval 
EQ.O - default set to 1 

0) 

6-10 

NPB 

Number of displacement print blocks 
EQ.O - print all nodal components 
GT.O - print according to table 
(NPB.LE.8) 

(2) 

11-15 

IDC 

Displacement printout code 
EQ.O - no displacement printout 
EQ.l - print displacements 

(3) 

16-20 

IVC 

Velocity print out code 
EQ.O - no velocity printout 
EQ.l - print velocities 

(3) 

21-25 

I AC 

Acceleration print out code 
EQ.O - no accelerations printed 
EQ.l - print accelerations 

(3) 

26-30 

IPGN 

Controls printing of generated nodal 
data 

EQ.O - print all generated nodes 
EQ.l - no print out 

(4) 

31-35 

IPGL 

Controls printing of generated nodal 
force components 

EQ.O - print all generated nodal force 
components 
EQ.l - no print out 

(5) 

36-40 

IRITRS 

Restart tape flag 

EQ.O - no restart tape written 

EQ.l - Restart tape written 

(6) 

41-45 

IPSM 

Controls printing of storage information 
messages 

EQ.O - print all messages 
EQ.l - no print out 

(7) 


NOTES 

(1) Program results will be printed every IPRI steps. Additional print 

options controlling output of displacements, velocities, accelerations, 
strains and stresses are specified by the variables NPB, IDC, IVC, and 


26 



II. PROGRAM CONTROL CARDS (continued) 


IAC on this card and on the element cards (Sections V.l or VI. 2). 

(2) Print blocks may be used to selectively print blocks or groups 
of nodal displacements (and velocities and accelerations). Each 
block of nodes is defined by the node number of the first and 
last node in the block (see card II. 5) 

(3) If NPB.EQ.O, then all displacements are printed regardless of the 
values of IDC, IVC and IAC. If NPB.NE.O, then IDC, IVC, and IAC 
must be input as 1 if displacements, velocities and accelerations, 
respectively, are to be printed. 

(4) When using the substructure input of nodal and element data 
(MGO = 1), this flag may be used to suppress the output of 
nodes automatically generated by the program. It is strongly 
suggested that the generated nodes be output and checked (for at 
least the first run of a problem). 

(5) When pressure loads are input, the program transforms these to 
generalized nodal forces. This flag may be used to suppress 
the output of these nodal forces. 

(6) A restart tape(s) which contains solution information for the final 
solution step may be written for future use in a restart of a 
problem. Information regarding these restart tapes and the setup 
of a restart is contained in Section PROGRAM RESTART. 

(7) The program uses "dynamic storage allocation" to allocate space for 
several key variables in common block AMAT. During the course of 

a solution, messages are printed which reflect how much storage 
is being used for each phase of the solution. 
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II. PROGRAM CONTROL CARDS (continued) 

II. 5 Print Block Table (16I5}*0MIT THIS CARD IF NPB.EQ.O 


col umns 

variable 

description 

notes 

1-5 

I PNODE ( 1 ,1 ) 

First node of print block no. 1 

(1) 

6-10 

IPN0DE(2,1 ) 

Last node of print block no.l 


11-15 

IPNODE ( 1 ,2) 

First node of print block no. 2 



etc. 



NOTES 




(1) The entries for each 
first and last node 

print block are input in pairs, i.e., 
of the print block. Nodal information 

the 

is 


printed at all nodes having node numbers between each pair of 
node numbers forming the print block. IF NPB.EQ.O, OMIT THIS 
CARD . 
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II. PROGRAM CONTROL CARDS (continued) 


II. 6 Mass And Damping Control Parameters* (41 5, 2F10.0)-0MIT THIS CARD 
IF IMASS. EQ.O 


col umns 

variables 

description 

notes 

1-5 

I DAMP 1 * 

Flag indicating damping type 
EQ.O - no Rayleigh damping 
EQ.l - Rayleigh damping 

(1) 

6-10 

IMASSN 1 ' 

Number of concentrated nodal 
masses 

(2) 

11-15 

IDAMPN f 

Number of concentrated nodal dampers 

(2) 

16-20 

IOPE 

Time integration method used 
EQ.O - default reset to 1 
EQ.l - Wilson's theta method 
EQ.2 - Newmark 1 s method 

(3) 

21-30 

OP VAR (1) 

First integration parameter 

1) If IOPE. EQ.l, parameter is e 
EQ.O - default reset to 1.4 

2) If IOPE. EQ.2, parameter is y 
EQ.O - default reset to 0.5 

(3) 

31-40 

0PVAR(2) 

Second integration parameter 

1) if IOPE. EQ.l, parameter is not 
applicable 

2) if IOPE. EQ.2, parameter is 6 
EQ.O - default reset is specified 
in note (3) 

(3) 


NOTES 

(1) In a dynamic analysis (IMASS. GT. 0) , Rayleigh damping may be included. 
The damping coefficients are defined in Section XII. 

(2) Concentrated nodal masses and dampers may be added if desired. These 
are input in Section XIV. 

(3) For a dynamic analysis, the time integration method must be selected. 
The variables OPVAR(l) and 0PVAR(2) contain the values of the para- 
meters required in the time integration method specified by IOPE. 

The Newmark e-method requires one parameter (e) which is normally set 
to 1.4. 

The Wilson 6-method uses the paremeters y and 6. For linear analyses 
y = 1/2 gives rise to no numerical damping and 6 = 1/4 corresponds 
to constant acceleration over a time step. For linear, unconditional 
stability, then y > 1/2 and p > (y+ l/2) 2 /4. If 6 is not specified 
(input as zero), then 6 is set to the default value given by the 
linear, stability criterion. 


f Not currently operational; must be input as zero. 
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II. PROGRAM CONTROL CARDS (continued) 


II. 7 Displacement Incrementation Variables (215, F10.0)-0MIT THIS CARD 
IF IDINC.EQ.O 


columns 

variable 

description 

notes 

1-5 

ININC 

Node number where degree-of-freedom 
is to be incremented 

0) 

6-10 

IDIR 

Di spl acement di recti on 
EQ.l - X direction 
EQ.2 - Y direction 
EQ.3 - Z direction 

(2) 

11-20 

DINC 

Displacement increment 

(3) 

NOTES 

(1) The global node 

number where a displacement increment is to 

sped fied 

/llAA 1 \ 


is required. If the substructure option for mesh generation (MG0=1) 
is being used, the required node number will not be known until a 
data check run is made and the automatically generated node numbers 
printed by the program. 

(2) IDIR refers to the global direction of the degree-of-freedom at a 
node. 

(3) The magnitude of the displacement increment is required; it may be 
positive or negative. 
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NODE AND ELEMENT INPUT 


The input data required to specify the node numbers, nodal coordinates, 
element numbers, element-node connectivity and element material properties 
is different depending upon which input option (MGO) is selected. 

If the individual nodal coordinates and element-node connectivity 
(MG0=0) are to be input, then the instructions in Sections III-VI should be 
followed. 

If the nodal and element geometry are being automatically generated 
using the substructure option (MG0=1 ) ,, then the input instructions in 
Section VII should be followed. 
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III. NODAL POINT DATA (MG0=0) 


The input instructions in Section 1 1 1 -V I are used only if MG0=O. 


III . 1 

Node and Element Group Data (315) 


column 

variable 

description 

notes 

1-5 

NUMNP 

Total number of nodal points 

(1) 

6-10 

NEGL 

Number of linear element groups 
EQ.O - all elements are nonlinear 

(2) 

11-15 

NEGNL 

Number of nonlinear element groups 
EQ.O - all elements are linear 

(2) 

NOTES 





(1) NUMNP is the total number of nodes to be defined in Section 
III. 2. 


(2) Elements are input by groups each of which have like characteristics 
(same element type, material model, type of geometric nonlinearity , 
and integration order). Linear elements are geometrically and 
materially linear and hence their' stiffness matrices are constant 
and formed only once. Nonlinear elements have their stiffness 
properties continually updated during the solution. 

The elements may be divided into groups for three reasons: (1) for 

input convenience - A particular group of elements may lend itself 
to being called a group simply because of the geometry of the structure, 
because of the node numbering, etc. Each element is assigned ascending 
element numbers from 1 to the total number of elements in that group. 

(2) There may be several regions each of which have different material 
properties, different element types,' etc. and each has to be treated 
as a different group. (3) Conservation of dynamic storage. Each 
element group uses repeatedly the same storage space in common block 
AMAT and hence less space is required if the the total number of 
elements is divided into several groups. 

For computational economy, all 'linear element groups are input together 
(first) and all nonlinear element groups are input together (after the 
linear element groups). 
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III. NODAL POINT DATA (continued) 


III. 2 Nodal Point Data (A1, 14, A1 , 14, 515, 3F10.0, 15) 


column 

1 


2-5 


variable description 

CT Symbol describing the coordinate 

system for this node 
EQ. - (blank) cartesian (X,Y,Z) 
EQ.X - X-cylindrical 

N Node (joint) number 

GE.l and LE.NUMNP 


notes 

0 ) 


( 2 ) 


6 PSF Print suppression flag (ignored (3) 

unless N.EQ.l ) 

EQ. - (blank) no suppression 
EQ.A - Suppress ordered list of node 
coordi nates 

EQ.B - suppress list of equation 
numbers 

EQ.C - both A and B, above 


7-10 I D( 1 ,N) 
11-15 ID(2,N) 
16-20 ID(3,N) 


X-translation boundary 
Y-translation boundary 
Z-translation boundary 


code (4) 

code 

code 


21-25 ID(4,N) X-rotation boundary code 
26-30 ID(5,N) Y-rotation boundary code 
31-35 ID(6,N) Z-rotation boundary code 


36-45 

X(N) 

X (or z)-coordinate 

46-55 

Y(N) 

Y (or r)-coordinate 

56-65 

Z(N) 

Z (or e)-coordinate (degrees) 


( 5) 


66-70 KN Node number increment for node (6) 

data generation 
EQ.O - no generation 


NOTES 

(1) A special cylindrical coordinate system is allowed for the speci- 
fication of node coordinates. If an "X" is entered in card column 
one (1), then the entries in cc 36-65 are with respect to a cylin- 
drical (r,8,z) system rather than the standard cartesian (X,Y,Z) 
system. The generator (i.e., the z-axis) of the "X-cylindrical" 
system coincides with the global X-axis, and e is zero (0) in the 
global X-Y plane (see figure below). 
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III. NODAL POINT DATA (continued) 


z-axis e-zero plane TRANSFORMATIONS 


X X-Y X = z 

Y = r cose 
Z = r sine 



x 


Cylindrical input is merely a user convenience for locating nodes 
in the (X,Y,Z) system, and no other references to the cylindrical 
system are implied; i.e., boundary conditions, output displacements 
etc. are still referenced to the (X,Y,Z) system. 

(2) Nodal data must be defined for all (NUMNP) nodes. Node data may 
be input directly (i.e., each node on its own individual card) or 
the generation option may be used if applicable (see note 6, below). 
Admissible node numbers range from "1" to the total number of nodes 
(NUMNP). Node numbers may not be repeated or omitted. The last 
node that is input must be "NUMNP". 

(3) The print suppression flag (PSF) is used to eliminate the second 
printing of ordered node coordinates or to suppress printing of 
equation number assignments (or both). The PSF character. is entered 
on the card for node one (1) only. 

(4) Boundary condition codes can only be assigned the following values 
(M=l ,2, . . . ,6) ~ 


M 



III. NODAL POINT DATA (continued) 


ID(M,N) = 0; unspecified (free) displacement 

ID(M,N) = 1; deleted (fixed) displacement 

ID(M,N) = -1; deleted (fixed) displacement used for 
generation 

An unspecified [ID(M,N)=0] degree of freedom is free to translate 
or rotate as the solution dictates. Concentrated forces (or 
moments) may be applied in this degree of freedom. 

One system equilibrium equation is required for each unspecified 
degree of freedom in the model. The maximum number of equilibrium 
equations is always less than six (6) times the total number of 
nodes in the system. 

Deleted [ID(M,N)=1] degrees of freedom are removed from the final 
set of equilibrium equations. Deleted degrees of freedom are used 
to define fixities (points of external reaction), and any loads 
applied in these degrees of freedom are ignored by the program. 

Nodes that are used for geometric reference only (i.e., nodes not 
assigned to elements) must have all six (6) degrees of freedom 
deleted. Nodal degrees of freedom having undefined stiffness (such 
as rotations in the TRUSS model, out-of-plane components in a two- 
dimensional model, etc.) must be deleted. 

Independent of the actual entries posted for the ID(M,N) in cc 7-35 
of the Nodal Data, any master boundary condition deletions (i.e., 

I DOF(M) . EQ . 1 , M-l,2,...,6) which are given in cc 1-6 of the card in 
Section II.l will be used for all nodes. Suppose that all rotations 
X,Y,Z have been deleted by means of the master codes (i.e., IDOF(M). 
EQ.l , M-4,5,6), then ID(MN ,) (M=4,5,6) will be set to "l", and data 
in cc 21-35 is ignored by the program. 

Two dimensional elements must use the Y-Z plane, i.e. plane strain, 
plane stress, and axi symmetric elements can only be defined in the 
Y-Z plane, and nodal degrees of freedom having undefined stiffness 
must be deleted. 

For generation purposes the value ID(M,N)=-1 can also be used. In 
this case, if the corresponding value on the next input card is zero 
(0), it is set equal to "-1 11 . Considering the deletion of degrees 
of freedom a minus one (-1) has the same meaning as a plus one (+1). 

(5) For the case CT (cc 1) equal to the character "X", the data input in 
cc 36-65 are interpreted as the cylindrical (r, 9, z) coordinates of 
node "N" . The table in notef 1) contains the formulae used by the 
program to compute the cartesian (X,Y,Z) coordinates of node N from 
the cylindrical coordinate values given in cc 36-65. The origin of 
cartesian and cylindrical systems is the same point. Cylindrical 
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III. NODAL POINT DATA (continued) 


coordinate data (if input on the card) are printed as read in the 
first list of node coordinates. The second listing of coordinates 
(if not suppressed with a PSF or "A" or "C") is an ordered printing 
of all nodes with coordinates converted to the (X,Y,Z) system. 

Note that boundary condition codes always reference the (X,Y,Z) 
system even if the node happens to be located using the cylindrical 
coordinate option. 

(6) Node cards need not be input in node order sequence; eventually, how- 
ever, all nodes in the set [1, NUMNP] must be defined. Node data 
for a series of nodes 

[N-j , N l +1*KN 1 , N 1 +2*KN 1 , ..., Ng] 

may be generated from information given on- two (2) cards in sequence -- 
CARD 1 - CT 1 ,N 1 ,ID(N 1 ,1),...,ID(N 1 ,6),X(N 1 ),...,KN 1 

CARD 2 - CT 2 ,N 2 ,ID(N 2 ,1),... ,ID(N 2 ,6),X(N 2 ),...,KN 2 

KN-, is the node generation parameter given on the first card in the 
sequence. The first generated node is N-j+1*KN-t; the second generated 
node is N-,+2*KNj , etc. Generation continues until node number 
N«-KN, is established. Note that the node difference N ? -N-. must be 
evenly divisible by KN-j. 

In the generation the boundary conditions codes (ID(L,J)values) of 
the generated nodes are set equal to those of node N-,. The coordi- 
nate values (cartesian (X,Y,Z) or cylindrical (r, e, z) are inter- 
polated linearly. Note that coordinate generation is only possible 
if both nodes N-jand N 2 are given in the same coordinate system. 

As an example of the use of the node generation parameter, consider 
the sketch below with nodes labeled as shown: 


5 

9 

13 

17 

21 


6 

10 

14 

18 

22 


7 

11 

15 

19 

23 


8 

12 

16 

20 

24 


25 

26 

27 

28 
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III. NODAL POINT DATA (continued) 


The nodal data may be generated using the following input: 


Card 1 
Card 2 
Card 3 
Card 4 
Card 5 
Card 6 
Card 7 


N=l, KN=4 


N = 25, KN = 0 


N = 2, KN = 4 
N = 26, KN = 0 
N = 3, KN = 4 
N = 27, KN = 0 
N = 4, KN = 4 


Generates nodes 5, 9, 13, 17, 21 
Generates' nodes 6, 10, 14, 18,' 22 


Generates nodes 7, 11, 15, 19, 23 


Card 8 


N = 28, KN = 0 


Generates nodes 8, 12, 16, 20, 24 


In using the above, one must be careful so that the generated boundary 
codes are the ones that are desired. In the above example, the 
boundary conditions specified on card 1 will be applied to the 
generated nodes, 5, 9, 13, 17, and 21. 
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IV. TRUSS ELEMENTS* 

TRUSS elements are two-node members allowed arbitrary orientation 
in the X,Y,Z system. The TRUSS transmits axial force only, and in 
general is a six (6) degree of freedom element (i.e., three global 
translation components at each end of the member, see Fig. 4). 


IV.l. 

Element Control 

Card (2014) 


col umn 

variable 

description 

notes 

1-4 

NPAR(l) 

Enter the number "1" 


5-8 

NPAR(2) 

Number of TRUSS elements in this group; 
GE.l 

(1) 

9-12 

NPAR(3) 

Type of geometric nonlinear analysis 
EQ.O - default set to "1" 

EQ. 1 - geometrically linear 
EQ.2 - geometrically nonlinear (updated 
formulation) 

(2) 

Lagrangian 

56-60 

NPAR(15) 

Material model number 
EQ.l - linear elastic 
EQ.2 - nonlinear elastic 

(3) 

61-64 

NPAR(16) 

Number of different sets of section/ 

material properties 

GE.l 

EQ.O - default set to"!" 

(4) 

65-68 

NPAR(1 7) 

Number of material model constants per 
set 

EQ.O if NPAR (15). EQ.l 
GE.4 if NPAR (15).’ EQ.2 

(5) 


NOTES 

(1) TRUSS element numbers begin- with one (1) and end with the total 
number of elements in this group, NPAR(2). Element data are input 
in Section IV. 4, below. 

(2) The paremeter NPAR{3) is applicable only if the element group is 
nonlinear. If NPAR(3) .EQ. 1 , no geometric nonlinearities are taken 
into account, i.e. the geometric stiffness matrix is not included. 

If NPAR(3).EQ.2 large displacement effects are included in the 
analysis, but small strains are assumed in the calculation of element 
forces. 

(3) In any one element group only one material model can be used, and 
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Fig. 4 Truss Model Parameters 
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IV. TRUSS ELEMENTS (continued) 


this model type is defined by the entry NPAR(15). If NPAR(15).EQ.l 
the model is defined by Young's modulus only and NPAR(17).EQ.O. If 
NPAR(15).EQ.2 the stress-strain curve is input in Section IV. 3 below. 

The model defined for the element group must be consistent with the 
nonlinear formulation used (defined by NPAR(3)) and the requirement 
of equilibrium iteration as defined on the Program Control Cards 
(Section II). As stated previously, equilibrium iterations can only 
be performed if the model allows for iteration, and if at least one 
nonlinear element group is used in the analysis. 


(4) The variable NPAR(16) defines the number of sets of material/section 
properties to be read in Sections IV. 2 and IV. 3 below. 

(5) NPAR(17) is the variable NCON used in Section IV. 3(c). 
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IV. TRUSS ELEMENTS (continued) 


IV. 2. Linear Elastic Material /Section Property Cards 

Skip this set of cards if NPAR(15).NE.l . 
Otherwise read NPAR(16) sets of cards. 


IV. 2. 

a. Material 

Number Card (15) 


column 

variable 

description 


1-5 

N 

Material /section number 
GE.li and LE.NPAR(16) 


IV. 2. 

b. Property 

card (4F10.0) 


col umn 

variable 

description 

notes 

1-10 

E(N) 

Young's modulus 

0) 

11-20 

AREA(N) 

Cross-sectional area 


21-30 

DEN(N) 

Mass density (for dynamic analysis) 


31-40 

STRAI(N) 

Initial strain (axial) 



NOTES 

(1) NPAR(16) different linear elastic materials are input in this section 
provided NPAR(15) .EQ. 1 . Note that one material /section is defined to 
have the same Young’s modulus, area, mass density and initial strain. 
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IV. TRUSS ELEMENTS (continued) 

IV. 3. Nonlinear Elastic Material /Section Property Cards 

Skip this set of cards if NPAR(15) .NE.2, otherwise read 
NPAR(16) sets of cards. 


IV. 3. a. 

Material number card (15) 


col umn 

variable 

description 


1-5 

N 

Material/section number 
GE.l and LE.NPAR(16) 


IV. 3. b. 

Section property card (3F10.0) 


column 

variable 

description 

notes 

1-10 

AREA(N) 

Cross-sectional area 

0) 

11-20 

DEN(N) 

Mass density for dynamic analysis 


21-30 

STRAI(N) 

Initial axial strain 


IV.3.C. 

Stress-strain 

curve card (8F10.0) 


column 

variable 

description 

notes 

1-10 

PROP ( 1 ,N) 

■j 

Strain at point 1, e 

(2) 

11-20 

PR0P(2,N) 

2 

Strain at point 2, e 



PROP (JnICON /2 s N)Strain at point NC0fV2 
PR0P(NC0N/2+l ,N) Stress at point 1, 
PROP(NCON/2+2 ,N) Stress at point 2, 2 


PR0P(NC0N,N) Stress at point NCON/2 
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IV. TRUSS ELEMENTS (continued) 
NOTES 


(1) One section property card is defined to have the same area, density 
and initial strain. 

(2) The stress-strain curve is defined by straight lines between the 
input points (e 1 , 0 1 ). From the stress-strain curve total stresses 
and the tangent modulus are evaluated for a given strain (see 
Figure 4). 

The variable NCON was defined in Section IV. 1 by the variable 
NPAR(17). 

This model can only be used in a nonlinear element group. 
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IV. TRUSS ELEMENTS (continued) 


IV. 4. Element Data Cards (415, F10.0, 315) 

NPAR(2) elements must be input and/or generated in this 
section in ascending sequence beginning with "1". 


col umn 

vari able 

des cri pti on 

notes 

1-5 

M 

TRUSS element number; 
GE.l and LE.NPAR(2) 


6-10 

II 

Node number at one end 

0) 

11-15 

JJ 

Node number at other end 
GE.l and LE.NUMNP 


16-20 

MTYP 

Material property set number; 
GE.l and LE.NPAR(16) , 


21-30 

PINIT 

Initial axial force in the TRUSS: 
GT.o, tension 

( a) 

31-35 

IPS 

Flag for printing axial stress in 
TRUSS element; 

EQ.O, no printing 
EQ.l, print element stress 


36-40 

KG 

Node generation increment used to 
compute node numbers for missing 
elements; 

EQ.O - default set to"l" 

(3) 


NOTES 

(1) Refer to Figure 4. 

(2) PINIT is the axial force in the TRUSS at zero node displacements, 
and zero initial strain in the truss. The initial strain defined 
for the material property set number gives rise to an additional 
force in the truss. 

(3) Elements must be input in increasing element number order. If cards 
for elements [M+l, M+2,...,M+J] are omitted, these "J" missing 
elements are generated using MTYP and PINIT of element ”M" and by 
incrementing the node numbers of successive elements with the value 
of "KG"; KG is taken from the first card of the element generation 
sequence (i.e., from the "M-th" element card). 
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V. 2/D CONTINUUM ELEMENTS 

2/D CONTINUUM elements are 4- to 8-node isoparemetri c quadrilaterals 
which must be input in the global Y-Z plane. Figure 5 shows some 
typical 2/D CONTINUUM elements, the element can represent either 
planar (plane stress or plane strain condition) or axisymmetric 
solids, as illustrated in Fig. 6 . In both cases, each element node 
has two (2) translational degrees of freedom. 

When the element is used to represent an axisymmetric solid or shell, 
the global Z-axis is the axis of revolution. All elements must be 
located in the +Y half-plane. These conventions are illustrated in 
Fig. 6. 
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AXISYMMETRIC FINITE ELEMENT MODEL OF A RING 



Fig. 6 Typical Two-Dimensional Applications* 


V. 2/D CONTINUUM ELEMENTS (continued) 


V.l. Element Group Control Card (2014)* 


columns variable 

1-4 NPAR(l) 

5-3 NPAR(2) 

9-12 NPAR(3) 

17-20 NPAR(5) 

25-28 NPAR( 7) 

37-40 NPAR(IO) 

49-52 NPAR(13) 

57-60 NPAR(15) 


description notes 

Enter the number "2“ 

Number of 2/D CONTINUUM elements in (1) 

this group 

GE.l 

Code indicating type of geometric nonlinear (2) 
analysis 

EQ.O - default set to "l” 

EQ.l - linear (small strain) 

EQ.2 - nonlinear (large strain), total Lagrangian 
Eq.3 - nonlinear (large strain), updated 
Lagrangian 

Element type code 
EQ.O - axi symmetric 
EQ.l - plane strain 
EQ.2 - plane stress 

Maximum number of nodes used to (3) 

describe any one element 

GE.4 and LE.8 

EQ.O - default set to "8" 

Numerical integration order to be (4) 

used in Gauss quadrature formulae 

GE.2 and LE.4 

EQ.O - default set to "2" 

Number of stress output location tables (6) 
EQ.O - print stresses at integration points 

Material model number (6) 

GE.l and LE.12 
EQ.l - linear isotropic 
EQ.2 - linear orthotropic 
EQ.3 - variable tangent moduli 
EQ.4 - curve description model 
EQ.5 - curve description model with tension 
cut-off (plane strain only) 

EQ.6 - elastic-plastic (von Mises with 

combined kinematic-isotropic hardening) 
EQ.7 - elastic-plastic (Drucker-Prager) 
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V. 2/D CONTINUUM ELEMENTS (continued) 


column variables 


61-64 NPAR(16) 
65-68 NPAR(1 7) 
69-72 NPAR (18) 


description notes 

EQ.8 - incompressible nonlinear elastic 
(Mooney-Ri vl in/plane stress only) 

EQ.9 - incompressible nonlinear elastic 
(Mooney-Ri vl in/plane strain or 
axi symmetric) 

EQ.10 - elastic-plastic (von Mises with 
mechanical sublayer hardening) 

GE.ll - user supplied model 

Number of different sets of material (6) 

properties 

GE.l 

Number of constants per property set (6) 

EQ.O - if NPAR(15).LE.10 

GE.O - if NPAR(15).G7.11 

Dimension of storage array required for (6) 

element history 

EQ.O - If NPAR(15) .LE. 10 

GE.O - If NPAR(15). GE.ll 


NOTES 

(1) 2D/C0NTINUUM element numbers begin with one (1) and end with the 
total number of elements in this group, NPAR(2). Element data are 
input in Section V.4 below. 

(2) NPAR(3) is applicable for nonlinear element groups only and de- 
termines if geometrical nonlinearities are to be included in the 
analysis. If NPAR(3).EQ.l displacements and strains are assumed 

to be infinitesimal. In the total Lagrangian and updated Lagrangian 
formulations all geometric effects (large displacements and large 
strains) are included in the analysis. 

(3) NPAR(7) limits the number of nodes that can be used to describe any 
of the elements in this group. A minimum of 4 and a maximum of 8 
nodes are used to describe the 2/D CONTINUUM elements. For example, 
for the set of elements shown in Fig.V.l, NPAR(7) would be set to 
"7". Constant strain triangles are obtained by having nodes 3 and 4 
coinci de. 

(4) For rectangular elements an integration order of "2" is sufficient. 

If the 'element is distorted, a higher integration order need be used. 
Notice that apart from the larger computational effort in the cal- 
culation of the element matrices, more working storage may be re- 
quired if a nonlinear material model is used (see note (6) below). 

The consistent mass matrix is always calculated with an integration 
order of 3. 
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V. 2/D CONTINUUM ELEMENTS (continued) 


(5) Element stresses are calculated at the points defined in the stress 
output location table assigned to the element (see SEction V. 3.). 
NPAR(13) defines the total number of stress output location tables 
input in Section V.3. NPAR(13).EQ.O for all nonlinear element 
models, i.e. MODEL. GT. 2 in Section V.2. 

(6) Only one material model (defined by the value of NAP R( 1 5) ) is al- 
lowed in an element group. IF NPAR(15) is 1,2,... ,10, the model 
exists in the current library for 2/D CONTINUUM elements, and the 
entries for NPAR(17) and NPAR( 1 8) are ignored by the program. If, 
however, NPAR(15) is “11“ or larger, material constants must be 
read into a storage array, "PROP", which has dimensions NPAR(17), 
NPAR( 16) , i.e., property constants are stoned as, 

((PR0P(I,<3),I=1 ,NPAR( 17) ),«!=! ,NPAR(16)) 

Non-library models require that NPAR(17) be given (input) as a 
positive integer. In addition, NPAR(18) must be specified for 
non-standard nonlinear models for use in allocating storage in a 
working array, "WA", which contains element history. The NPAR(18) 
parameters which characterize element history must be retained for 
materially nonlinear elements so that properties can be chosen in 
the current solution step. The working storage array (WA) is dim- 
ensioned IDWA by NPAR(2); where 

IDWA * NPAR(18)*(NPAR{10)**2) 

NPAR(2) = number of elements 

i.e. for each integration point of each element, NPAR(18) storage 
locations are preserved during the solution. 

The model defined for the element group must be consistent with the 
nonlinear formulation used (defined by NPAR(3)), the element type 
(defined by NPAR(5}), and the requirement of equilibrium iteration 
as defined in Section II. As stated previously, equilibrium it- 
erations can only be performed if the model allows for iteration, and 
if at least one nonlinear element group is used in the analysis. 
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V. 2/D CONTINUUM ELEMENTS (continued) 

V.2. Material Property Data 

NPAR (16) sets of cards must be input in this section. Card "a" 
(material number card) is the same for all material models, but 
card(s) "b" ("material property card(s)' 1 ) depend on the material 
model number (NPAR(15)). 


V.2. a. 

Material number card (I5,F10.0) 


columns 

variable 

descripti on 

notes 

1-5 

N 

Material property set number; 
GE.l and LE.NPAR(16) 


6-15 

DEN(N) 

Mass Density of the material used in 
calculation of the mass matrix, p; 
GE.0.0 

(1) 

NOTES 





(1) The mass density defined is used directly in the calculation of 
the element matrix, i.e., no acceleration constants are applied 
to the variable DEN(N). 


V.2.b. Material property card(s) (8F10.0) 

MODEL 1 - linear isotropic* 

columns variable description 

1^10 'PROP (1 ,N) Young's modulus, E 

11-20 PR0P(2,N) Poisson's ration v 


NOTES 

(1) MODEL 1 is a linear library model defined by two (2) positive constants 
(E»v) ; i.e., if NPAR{15).EQ.1 , NPAR(17) is set to "2" by default. 

MODEL 1 can be used with linear or nonlinear element groups. Since 
the material constants are independent of history, NPAR{18) is set to 
"O' 1 by default. 


notes 

( 1 ) 
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X. 2/D CONTINUUM ELEMENTS (continued) 
MODEL 2 - linear orthotropic 


col iimn 

vari able 

description 

notes 

1-10 

PR0P(1,N) 

a-di recti on modulus, E 

a 

( 1 ) 

11-20 

PR0P(2,N) 

b-di recti on modulus, E^ 


21-30 

PRQP(3,N) 

c-di recti on modulus, E c 


31-40 

PR0P(4,N) 

strain ratio, u , 


41-50 

PR0P(5,N) 

strain ratio, v, 
ac 


51-60 

PR0P(6,N) 

strain ratio, v u 
be 


61-70 

PR0P(7,N) 

shear modulus, G ab 



NOTES 


MODEL 2 is linear (i.e., NPAR(18).EQ.O by default). The model 
requires seven (7) property constants; i.e., NPAR{17).EQ.7 by 
default if NPAR(15) .EQ.2. The constants(E., E, ,...,G h ) are defined 
in a meterial coordinate system. (a, b,c) in a which "c" looresponds 
to the out-of-plane direction (or global X), and the a-b-plane 
coincides with the plane of the elements (Y-Z) . a,b are principal 
material axes, and 3 is the angle (in degrees) between element 
edge 1-2 and the material +a-axis (Figure 7 ). The material 
angle 6 is read individually for each element in Section V.4, 
below. MODEL 2 can be used with linear or geometrically nonlinear 
elements. The moduli E a , Eb, E c , and G a b must all be greater than 
zero . 
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Fig. 7 Principal In-Plane Material Axes Orientation 
for Linear Orthotropic Material Model 2* 




X. 2/D CONTINUUM ELEMENTS (continued) 
MODEL 3 - variable tangent moduli* 


column 

1-10 

11-20 

21-30 

31-40 

41-50 

51-60 

61-70 


variable 


description 

notes 

PROP ( 1 ,N) 

Unloading bulk modulus, K uN 

(1) 

PR0P(2,N) 

K 0 



PROP ( 3 S N) 

K 1 

Constants to define the loading 
bulk modulus 


PR0P(4,N) 

K 2 



PROP (5, N) 

G 0 



PR0P(6,N) 

T-j ( <0 ) 

Constants to define the loading 


PR0P(7,N) 

T-j (<0) 

and unloading shear moduli 



NOTES 

(1) MODEL 3 is a nonlinear material model requiring seven (7) constants 
(i.e., NPAR( 1 7) - Ey . 7 by default) per property set. Ten variables, 
namely thepreviously calculated stresses, strains, the maximum 
pressure ever reached and a loading parameter, need be stored per 
integration point, i.e., NPAR(18) = 10 for this model. 

The material model cannot be used in plane stress analysis. 


NOTATION 

For the notation used in the description of the nonlinear material _ 
models refer to the report "Static and Dynamic Geometric and Material 
Nonlinear Analysis," SESM Report No. 74-4, February 1974 (Ref, 3). 
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V. 2/D CONTINUUM ELEMENTS (continued) 

MODEL 4 or 5 - curve description model* 

a. Material moduli cards (supply NPAR(16) sets of cards) 


col umn 

varaible 


description notes 

first card 




1-10 

PR0P(1,N) 


volumn strain at point "l" (1) 

11-20 

PR0P(2jN) 

e v 

• • • • 

21-30 

PROP (3, N) 

e v 


31-40 

PR0P(4 5 N) 

< 


41-50 

PR0P(5,N) 

4 


51-60 

PROP (6, N) 

4 


61-70 

PR0P(7 5 N) 

K 1 

k ld 

loading bulk modulus at point "1" 

71-80 

PROP ( 8 a N) 

K 2 

k ld 

• » • • 

second card 



1-10 

PROP (9 ,N) 

K 3 

k ld 


11-20 

PR0P(10,N) 

k l!d 


21-30 

PR0P(11,N) 

k ld 


S' 1 -40 

PR0P(12 ,N) 

K 6 

k ld 


41-50 

PR0P(13,N) 

K 1 

k un 

unloading bulk modulus at point "1" 

51-60 

PROP ( 1 4 ,N ) 

i/2 

*UN 


61-70 

PROP (15 ,N) 

K 3 

k un 


71-80 

PROP(16,N) 

K 4 

k un 


third card 




1-10 

PR0P(1 7,N) 

K 5 

uN 

(2) 

11-20 

PR0P(18,N) 

k un 


21-30 

PROP ( 1 9 ,N ) 

<fn 

loading shear modulus at point "1" 


55 



V. 2/D CONTINUUM ELEMENTS (continued) 


col umn 

vari able 


31-40 

PROP(20,N) 

r 2 

g ld 

41-50 

PROP ( 2T ,N) 

r 3 

g ld 

51-60 

PROP (22,N) 

r 4 

g ld 

61-70 

71-80 

PR0P(23,N) 

PR0P(24,N) 

g ld 6 

4 D 


description notes 


NOTES 

(1) MODELS 4 and 5 are nonlinear material models for which NPAR(17) .EQ.24, 
by default. To characterize history, twelve (12) variables, namely, 
the previously calculated stresses and strains, the maximum volumetric 
strain ever reached in loading, ground pressure, the volume strain 
due to ground pressure, and the angle of "cracking" need be saved for 
each integration point, e.i. NPAR(18) . EQ. 12- 

Note that 


e J > 
V 


*v 1 5 K UN - K LD and G LD 


< i.5 k; 


LD 


where "J" is the J-th table input point. All six (6) input points must 
be defined, as shown in Fig. 8 . Linear interpolation is used to 
define the loading and unloading bulk moduli, ^K.r,, and the loading 

shear modulus, t G. n , at time t. Note that these three cards must be 
supplied NPAR(16) L times. 


(2) The unloading shear modulus is calculated as ( K UN , K, D 


) 4 G 


'i n* 


b. Tension cut-off definition card 

This card defines the parameters to be used in the tension cut-off model. 


i.e. , 

MODEL 5. 

Omit this card if MODEL. EQ. 4. 


column 

variable 

description 

notes 

1-10 

GAMMA 

material density 

(1) 

11-20 

21-30 

STIFAC 
SHE FAC 

normal stiffness reduction factor 
shear stiffness reduction factor 

(2) 
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V. 2/D CONTINUUM ELEMENTS (continued) 
NOTES 


(1) GAMMA is used to calculate the in-situ gravity pressure at the 
element integration points. It is assumed that the material 
provides a hydrostatic pressure distribution, calculated at an 
element integration point using 

N 

P = 2 h, p. 
i=l 1 1 

where h^ are the element interpolation functions, p. is the pressure 
at the element nodes, and N is the total number of element nodes. 

The pressure p-j is calculated using p-j = GAMMA * Z-j , where Zj is the co- 
ordinate of node i. This convention defines the assumed orientation 
of the Y-Z coordinate system. (All two-dimensional analyses must be 
carried out in the Y-Z plane). From the hydrostatic pressures and 
the bulk moduli input above, the in-situ volume strain at every inte- 
gration point is calculated. When finding bulk and shear moduli 
from the curve description shown in Fig. 8 , in-situ volume strain 
is included. That is, for MODEL 5 


t _ t P 
V ~ e V 


+ 



where ^ey is total volume strain (abscissa of curves in Fig. 8 ), 

t P t l 

e^ is the in-situ volume strain, and e^ is the volume strain 

due to applied loads. Note that for MODEL 4 ^e,. is assumed to equal 

4 - 1 * 


(L) MODEL 5 contains the assumption that the material cannot withstand 

tensile stresses. When the tensile principal stress at an integration 
point exceeds the in-situ gravity pressure, a crack forms perpendicular 
to this principal stress and the stiffness across the crack is reduced 
by the factor STIFAC. Also, the shear stiffness at that integration 
point is reduced by the factor SHEFAC, which is meant to model the 
loss in shear strength due to cracking. Model 5 is to be used 
with plane strain elements only and unloading is not permitted. 
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V. 2/D CONTINUUM ELEMENTS (continued) 

MODEL 6 - Combined Isotropic-Kinematic Hardening, von Mises yield condition 


column 

CARD 1 1-10 

11-20 
21-30 
31-40 

41-50 

51-60 

61-70 

71-80 


variable 

description 

notes 

PROP ( 1 ,N) 

Young's modulus, E 

O) 

PR0P(2 ,N) 

Poisson's ratio, v 

O) 

PROP (3,N) 

dumber of points in piecewise linear 
Stress Strain Curve, RNSSC 

(2) 

PR0P(4,N) 

Number of points in piecewise yield 
surface SIZE vs equivalent plastic 
curve, RNSIZE 

(3) 

PROP (5 »N ) 

Type of hardening rule being used, 
TP LAS 

(4) 

PROP (6, N) 

Ratio of isotropic hardening to 
kinematic hardening, RITK 

(5) 

PR0P(7,N) 

■ Value of stress at first point in stress 
strain curve, ol 

Value of total strain at first point in 
stress strain curve, &j 

(6) 

PR0P(8,N) 

■(6) 


CARDS 

2,3,4 1-10 


11-20 


PR0P(9,N) Value of stress at second point in stress 
strain curve, o 2 

PROP ( 1 0 ,N ) Value of strain at second point in stress 

strain curve, 


PROP (5+2* 
RNSSC,N) 


PROP (6+2* 
RNSSC,N) 


PORP ( 7+2* 
RNSSC,N) 


PROP (8+2* 
RNSSC,N) 


Value of stress at last point in stress 
strain curve, a RNSSC 


Value of strain at last point in stress 
strain curve, t rns$c 


Value of yield surface size at first point (7) 
in yield surface size vs equivalent plastic 
strain curve, 


Value of equivalent plastic strain at first 
point in yield surface size vs equivalent 
plastic strain curve, 
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V. 2/D CONTINUUM ELEMENTS (continued) 


For MODEL "6" (continued) 

column variable description notes 


PROP (5+2* 

RNSSC+2*RNSIZE 

,N<) Value of yield surface size at last 

point in yield surface size vs 
equivalent plastic strain curve* k resize 

PROP (6+2* 

RNSSC+2*RNSIZE 

,N) Value of equivalent plastic strain at 

last point of yield surface size vs 
equivalent plastic strain curve* ^r^sjze 

NOTES 

(1) MODEL 6 is nonlinear library material model for which NPAR(17) .EQ.26 
by default. NPAR(18) .EQ.15, since 15 history dependent variables 
must be stored for each integration point. 

The Young's modulus and Poisson's ratio define the initial elastic 
behavior of the material. 

(2) RNSSC is the number of points in a piecewise linear stress strain 
curve, and includes the end point of the last linear segment, but 
not the point STRESS=0, STRAIN=0. Maximum allowable value is RNSSC= 
5, the minimum value is RNSSC=2. 

(3) RNSIZE is the number of points in a piecewise linear yield surface 
size vs equivalent plastic strain curve, and includes both the end 
point of the last linear segment as well as the point STRESS=yield 
stress, EQPLAS=0. RNSIZE=0 for TPLAS.NE.3. 

(4) TPLAS EQ.O Isotropic hardening 

EQ.l Kinematic hardening (Prager-Ziegler) 

EQ.2 Combined isotropic, kinematic 

hardening, constant ratio 

EQ.3 Combined isotropic, kinematic 

hardening, variable ratio 

(5) RITK is the ratio of isotropic to kinematic hardening for TPLAS.EQ.2, 
i.e., with RITK. EQ.l being pure isotropic, RITK. EQ.O pure kinematic. 
If TPLAS. NE. 2, RITK is ignored. 
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V. 2/D CONTINUUM ELEMENTS (continued) 
For MODEL "6 1 ' (continued) 


NOTES 

(6) Values of STRESS and STRAIN are input in pairs for the piecewise 
linear stress strain curve using the format 8F10.0. A maximum 

of five (5) pairs of data is allowed. Input RNSSC pairs of data. 

The input stress-strain curve must be monotonicly increasing such 
that PR0P( 5+2*1 ).LE. PR0P(5+2*(I+1) ) and PROP (6+2*1) .LT. PROP 

( 6+2*( 1+1 ) ) for any I. See also Fig. 9. 

(7) Values of YSSIZE (Yield Surface SIZE) and EQPLAS (Equivalent PLAStic 
strain) are input in pairs for the piecewise linear yield surface 
size vs equivalent plastic strain curve. Data for this curve need 
not be input unless TPLAS.EQ.3. For TPLAS.LT.3, a compatable curve 
is generated internally in the program, appropriate for the hardening 
rule indicated by TPLAS. See also Fig. 9. 

NOTE: Regardless of the value of RNSSC and RNSIZE, four cards will 

be read on 8F10.0 format to generate the PROP array. Blank entries 
for a property value will be read and ignored after the last non- 
zero entry. FOUR CARDS MUST BE SUPPLIED. 
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Fig. 9 Stress and Yield Surface Size vs. Strain Curves 
for Combined Isotropic-Kinematic Hardening Model 


Grnsize € 
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V. 2/D CONTINUUM ELEMENTS (continued) 

MODEL 7 - elastic-plastic material, Drucker-Prager yield condition * 


column 

variable 

description 

notes 

1-10 

PROP ( 1 ,N) 

Young's modulus, E 

0) 

11-20 

PR0P(2,N) 

Poisson's Ratio, v 


21-30 

PR0P(3,N) 

Angle of friction, e 


31-40 

PR0P(4,N) 

Cohesion, C 


NOTES 





(1) MODEL 7 is a nonlinear library material model requiring four (4) 
constants (i.e., NPAR(17) .EQ.4 by default) per material property 
set. Eleven (11) parameters need to be stored per integration 
point, i.e. NPAR(18) .EQ.ll by default. The model assumes elastic 
perfectly plastic conditions. For a further description of the 
material model parameters, refer to the report "Static and Dynamic 
Geometric and Material Nonlinear Analysis," SESM Report No. 74-4, 
Civil Engineering Department, University of California, Berkely, 
California, February 1974 (Ref. 3). 
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V. 2/D CONTINUUM ELEMENTS (continued) 


MODEL 8 or 9 - incompressible 

nonlinear elastic material, Mooney-Ri vlin model 

col umn 

variable 

description 

notes 

1-10 

PROP(l.N) C 1 

C-, and C 2 are material constants 
in the Mooney-Ri vlin expression 

0) 

11-20 

PROP ( 1 ,N) C 2 

for strain energy. 


NOTES 




(1) MODEL 

8 and 9 are nonlinear library material models for which NPAR(17) 


is set to "2". Since the material is elastic, no history-dependent 
variables need be stored at integration points and NPAR(18) is set to 
"0". MODEL 8 is for plane stress only and MODEL 9 is for plane strain 
or axi symmetric only. C, and C 2 are material constants in the expression 
for strain energy: 

v = C ] (I 1 -3) + C 2 (I 2 -3) 

where I-| and 12 are the first and second strain invariants, respectively. 
A more complete description of the model is given in Ref. 5. 
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V. 2/D CONTINUUM ELEMENTS (continued) 



MODEL 10 

- Mechanical 

Sublayer Model , von Mises yield condition 



col umn 

vari able 

description 

notes 

CARD 1 

1-10 

PROP ( 1 ,N) 

Young's modulus, E 

(1) 


11-20 

PR0P(2,N) 

Poisson's ratio, v 

O) 


21-30 

PR0P(3,N) 

Number of points in piecewise linear 
Stress Strain Curve, RNSSC 

(2) 


' 31-40 

PR0P(4,N) 

Value of stress at first point in stress 
strain curve, 

(2) 


41-50 

PR0P(5,N) 

Value of strain at first point in stress 
strain curve, e-j 

(2) 



PR0P(8,N) 

Value of stress at third point in stress 
strain curve, a 2 


CARD 2 


PR0P(9,N) 

Value of strain at third point in stress 



strain curve, e 2 


PR0P(2+2*NSSC 

»N) Value of stress at last point in stress 

strain curve, c RNSSC 

PR0P(3+2*NSSC 

»N) Value of Strain at last point in stress 

strain curve , f_ k1cc _ 

RNSSC 

NOTES 

(1) MODEL 10 is a nonlinear library material model for which NPAR(-17) 

EQ.11 by default. NPAR(18) .EQ.20, since 20 history dependent 
variables must be stored for each integration point. 

The Young's modulus and Poisson's ratio define the initial elastic 
behavior of the material. 

(2) RNSSC is the number of points in a piecewise linear stress strain 
curve, and includes the end point of the last linear segment (unless 
that segment is perfectly plastic), but not the point STRESS=0, 
STRAIN=0. The maximum allowable value is RNSSC=4.0. If RNSSC. LE. 2.0, 
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v. 2/d CONTINUUM ELEMENTS (continued) 


a blank card B must still be supplied. It is also required that the 
input stress strain curve be monotonicly increasing so that PROP 
( 2+2*1 }.LT.PR0P(2+2*(I+1) and PR0P(3+2*I).LT.PR0P(3+2(I+1)). Note 
that the inequality reads .LT., not .LE. Figure 10 depicts a typical 
pi ecewise-1 inear stress-strain ctrrve. 



Fig. 10 Equivalent Uniaxial 5tress -Strain Curve for 
Mechanical Sublayer Model 
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V. 2/D CONTINUUM ELEMENTS (continued) 


MODEL II - user-supplied material model* 


column 

variable 

description 

notes 

1-10 

PR0P(1,N) 


(1) 

11-20 

PR0P(2,N) 



21-30 

PR0P(3,N) 



• 

• 

etc. 



NOTES 





(1) MODEL 11 is a non-standard model (i.e., is not in the material model 
library) and requires that property data be read 


(PR0P(J,N),J=1,NPAR(17),N=1 ,NPAR(16) 

where NPAR( 17) is user-supplied (i.e., no default value exists). 
Subroutines which allocate storage locations (in "WA") for element 
history (if the model is nonlinear) and which calculate the incre- 
mental stress-strain law and current stresses must be supplied. 

For a discussion on how to incorporate a user-supplied material 
model, see Ref. 2. 
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V. 2/D CONTINUUM ELEMENTS (continued) 

V.3. Stress Output Table Cards (915) * 

Skip this section of stresses at integration points are to be 
printed, i.e. NPAR(13) .EQ.O; otherwise supply NPAR(13) cards. 


column 

variable 



description 

1-5 

ITABLE(N,1) 

Stress 

output 

location 1 

6-10 

ITABLE(N,2) 

Stress 

output 

location 2 

11-15 

ITABLE(N,3) 

Stress 

output 

location 3 

16-20 

ITABLE(N,4) 

Stress 

output 

location 4 

21-25 

ITABLE(N,5) 

Stress 

output 

location 5 

26-30 

ITABLE(N,6) 

Stress 

output 

location 6 

31-35 

ITABLE(N,7) 

Stress 

output 

location 7 

36-40 

ITABLE(N,8) 

Stress 

output 

location 8 

41-45 

ITABLE(N,9) 

Stress 

output 

location 9 


notes 

( 1 ) 


NOTES 

(1) Stress tables are defined to provide flexibility in element stress 

output request's. Each element can refer to one of the tables defined, 
and the element stresses are then calculated at the points specified 
in the table. Refer to Fig. 11 for selection of the stress cal- 
culation points. The first "0" entry in a table will terminate that 
table. For example, if ITABLE(N,1) .EQ.7 and ITABLE(N,2) .EQ.O, then 
stresses will be printed at point 7 only whenever this stress table 
is referred to. 

The stress tables are only used when MODEL "1" or "2" define the 
material behavior. 
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STRESS PRINTOUT TABLE POINTS FOR 
2/D CONTINUUM ELEMENTS 


S 



3- POINT INTEGRATION 

PRINTING SEQUENCE FOR INTEGRATION 

POINT STRESSES FOR 2/D CONTINUUM ELEMENTS 


Fig. 11 Stress Print-Out Locations for Two-Dimensional Elements* 
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V. 2/D CONTINUUM ELEMENTS (continued) 

V.4. Element Data Cards (15, 13, 12, 2F10, 1015) * 

column variable description notes 

I- 5 M 2/D CONTINUUM element number; (1) 

GE.l and LE.NPAR(2) 

6-8 IEL Number of nodes used to describe this (2) 

element; 

EQ.O - default set to "NPAR(7)" 

LE.NPAR(7) 

9-10 IPS Number of the stress output table to be used (3) 

for stress calculations for this element (if 
a linear material model is used) or flag for 
printing stresses at element integration 
points (nonlinear material models; 

. EQ.O - no stress output for this element 

II- 20 BET Material angle, e, in degrees to be used 

in connection with material model "2" 

(linear orthotropic) 

21-30 THIC Element thickness; (4) ' 

GT.O - plane stress only 

31-35 MTYP Material property set number assigned 

to this element; 

GE.l and LE.NPAR(16) 

36-40 KG Node generation parameter used to compute (1) 

node numbers for missing elements 
(given on first card of a sequence); 

EQ.O - default set to "1" 

41-45 N0D(1) Global node number of element nodal point 1 (2) 

46-50 N0D(2) Global node number of element nodal point 2 

51-55 N0D(3) Global node number of element nodal point 3 

56-60 N0D(4) Global node number of element nodal point 4 

61-65 N0D(5) Global node number of element nodal point 5 

66-70 N0D(6) Global node number of element nodal point 6 

71-75 N0D(7) Global node number of element nodal point 7 

76-80 N0D(8) Global node number of element nodal point 8 
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V. 2/D CONTINUUM ELEMENTS (continued) 


NOTES 

(1) Elements must be input in ascending element number order. If data 
cards for elements [M+l , M+2,. . . ,M+J] are omitted, these "J" missing 
elements are generated using IEL, IPS, BET, THIC, and MTYP given on 
the card for element "M" and by incrementing node numbers of suc- 
cessive elements with the value /‘KG' 1 ; the value of KG used for incre- 
mentation is taken from the M-th element card, and only the non-zero 
nodes appearing on the M-th element card are incremented when gen- 
erating missing element data. * The first and last element must always 
be input. 


(2) The number of nodes in element "M" is defined by 11 lEL" . However, 
all 8 entries for NOD(I) are read from the element data card; if 
IEL.LT.8 the particular node locations not used in this element 
need be input as "0" in NOD(I). Fig. 12 defines the input 
sequence that must be observed for element node input. 

For example, the 6-node element (IEL.EQ.6) shown below 



is defined by 

NOD(I) =[XXXX0X0X] 
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V. 2/D CONTINUUM ELEMENTS (continued) 

NOTES 

(2 continued) 

where the nonzero entries (X) are the global mesh node numbers 
(Section III. 2) of the 6 nodes. 

(3) If IPS.EQ.O no stress output will be provided for this element. 

A nonzero entry for the element stress printout flag "IPS" has 

a different meaning when used with nonlinear material models than 
it does when used with linear material models: 

a) For linear material models (MODELS "1" and "2") IPS specifies 
the stress output table (Section V.3) to be used for this 
element. If no stress output tables were input (i.e., if 
NPAR(13).EQ.O) then IPS.EQ.1 specifies that stresses are to 

be printed at all integration points for this element, as shown 
in Fig. 11. 

b) For nonlinear material models ( MODELS. GT." 2") IPS.EQ.1 specifies 
that stresses are to be printed at all integration points for 
this element, as shown in Fig. 11. 

(4) The thickness of the element need only be defined in plane stress 
analysis. In plane strain analysis thickness is assumed to be "1". 
In axi symmetric analysis thickness is assumed to be one (1) radian. 


i-x 



VI. 3/D CONTINUUM ELEMENTS 

3/D CONTINUUM elements are 8- to 21-node isoparametric or sub- 
parametric curvilinear hexahedra. Typical elements are illustrated 
in Figure 13. 

These elements may be used to model three-dimensional solids or 
thick shells, two examples of which are shown in Figure 14;. 
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3/D CONTINUUM ELEMENT MODEL OF ARCH DAM 




VI. 3/D CONTINUUM ELEMENTS (continued) 


VI. 1. 

Element Group Control Card (2014)* 


col umn 

vari able 

descripti on 

notes 

1-4 

NPAR(l) 

Enter the number "3" 


5-8 

NPAR(2) 

Number of 3/D CONTINUUM elements in 

this group 

GE.l 

(1) 

9-10 

NPAR(3) 

Flag indicating type of geometric nonlinear (2) 
analysis 

EQ.O - default set to "1" 

EQ.l - linear (small strain) 

EQ.2 - nonlinear (large strain), total Lagrangian 
EQ.3 - nonlinear (large strain), updated 
Lagrangian 

25-28 

NPAR(7) 

Maximum number of nodes used to 
describe any one element 
EQ.O - default set to "21" 

(3) 

37-40 

NPAR(IO) 

Numerical integration order to be 
used in Gauss quadrature formula 
for r-s plane 
EQ.O - default set to "2" 

GE.2 and LE.4 

(4) 

41-44 

NPAR(ll) 

Numerical integration order to be 
used in Gauss quadratrue formula 
for t-di recti on 
EQ.O - default set to "2" 

GE.2 and LE.4 

(4) 

49-52 

NPAR(13) 

Number of stress output location 
tables 

EQ.O - print stresses at integration 
points 

(5) 

57-60 

■ NPAR( 15) 

Material model number 

(6) 


GE.l and LE.8 

EQ.l - linear elastic isotropic 
EQ.2 - linear elastic orthotropic 
EQ.3 - curve description model 
EQ.4 - elastic-plastic, combined kinematic 
isotropic hardening, von-Mises yield 
condi ti on 

EQ.5 - elastic-plastic, mechanical sublayer 
hardening, von-Mises yield condition 
EQ.6 - incompressible nonlinear elastic, 
Mooney- Ri vl 1 n mate ri a 1 
GE.7 - user supplied model 
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VI. 3/D CONTINUUM ELEMENTS (continued) 


col umn 

variable - 

description 

notes 

61-64 

NPAR(16) 

Number of different sets of material 

properties 

GE.l 

(6) 

65-68 

NPAR( 1 7) 

Number of constants per material 

property set 

EQ.O - if NPAR(15) .LE.@ 

GE.l - if NPAR(15.SE. 7 

(6) 

69-72 

NPAR(18) 

Dimension of storage array required 
for element history 
EQ.O - if NPAR( 15) .LE.6 
GE.O - if NPAR(15).GE.7 

(6) 


NOTES 

(1) 3/D CONTINUUM element numbers begin with one (1) and end with the 
total number of elements in this group, NPAR(2). Data for the in- 
dividual elements are input in Section VI. 4 below. 

(2) NPAR(3) is applicable for nonlinear element groups only and determines 
the type of geometrical nonlinearities to be considered in the 
analysis. 

(3) NPAR(7) is the maximum number of nodes that can be used to describe 
any one of the elements in this group, i.e. elements in this group 
must have less than or equal to NPAR(7) nodes. A minimum of 8 and 
a maximum of 21 nodes are used to describe 3/D CONTINUUM elements 

as indicated in Fig. 16 . Constant strain tetrahedra can be formed 
from 8-node elements by having nodes 1,2,3 and 4 coincide and nodes 
7 and 8 coincide. 

(4) The selection of appropriate integration orders depends on the 
element shape and strain state being considered. When the quantities 
being integrated vary irregularly a higher order is needed. An 
integration order of ll 2" is sufficient for most problems (especially 
those modeled with 16 node elements); however, an integration order Of 

3 should be used for all 20 node elements. Shell or plate problems often 
use thin elements in which strain varies more or less linearly through the 
thickness, but more irregularly in the plane of the surface of the shell. 
In these cases it is advantageous to specify a lower order for integration 
through the thickness, i.e., NPAR(ll) < NPAR(IO). 

The expense of stiffness formation is dependent on the integration order 
i.e.. 


n = NPAR(ll) * NPAR(IO) * NPAR(IO) 
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VI. 3/D CONTINUUM ELEMENTS (continued) 
NOTES 


where n is the number of points at which B T DB must be calculated 
in order to find the element stiffness by Gauss integration. 

The consistent mass matrix is always calculated using an integration 
order of three in each direction. 

(5) Element stresses are calculated at the points defined in the stress 
output location table (see SEction VI. 3) assigned to an element. 

NPAR(13) defines the total number of stress output location tables 
input in Section VI. 3. 

For nonlinear material models (NPAR(15).GE.3) stress output tables 
cannot be used, and stresses are printed at integration points; 

NPAR(13) is set to "0" in these cases. 

(6) Only one material model (defined by the value of NPAR(15)) is allowed 
in an element group. If NPAR(15) is 1, 3, 4, 5, or 6 the model exists 
in the current library for 3/D CONTINUUM elements, and the entries 
for NPAR(17) and NPAR(18) are ignored by the program. If, however, 
NPAR(15) is "7" or larger, material constants must be read into a 
storage array, “PROP", which has dimensions NPAR(17) by NPAR(16).; i.e., 
property constants are stored as, 

( ( PROP ( I ,J) , 1=1 ,NPAR(17)),J=1 ,NPAR(16)) 

Non-library models require that NPAR(17) be given (input) as a positive 
integer. In addition, NPAR(18) must be specified for non-standard non- 
linear models for use in allocating storage in a working array, "WA" 
which contains element history. The NPAR(18) parameters which charact- 
erize element history must be retained for materially nonlinear elements 
so that properties can be chosen in the current solution step. The 
working storage array(WA) is dimensioned IDWA by NPAR(2), where 

IDWA = NPAR( 18) * (NPAR(1G)**2) * NPAR(ll) 

NPAR(2) = number of elements 

i.e., for each integration point of each element NPAR(18) storage 
locations are preserved during the solution. 


The model defined for the element group must be consistent with the 
nonlinear formulation used (NAPR(3)) and the requirement of equilibrium 
iteration. Equilibrium iteration can be performed only if the material 
model allows for it and if there is at least one nonlinear element 
group. 
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VI. 3/D CONTINUUM ELEMENTS (continued) 

VI. 2. Material Property Data 

NPAR(16) sets of cards must be input in this section. Card "a" 
("material number card") is the same for all material models, but 
card(s) "b" ("material property card(s)") depend on the material 
model number NPAR(15). 


VI. 2. a. 

Material number card (I5,F10.0)* 


col umn 

variable 

description 

notes 

1-5 

N 

Material property set number 
GE.l and LE.NPAR(16) 


6-15 

DEN (N) 

Mass density of the material used in the 
calculations of the mass matrix 
GE.0.0 

(1) 


NOTES 

(1) The mass density defined is used directly in the calculationof the 
element mass matrix, i.e. no acceleration constants are applied to 
the variable DEN(N). 


VI. 2. b. Material property card(s) (8F10.0) 

MODEL 1 - linear isotropic elastic 

column variable description 

I- 10 P ROP ( 1 ,N) Young's modulus, E 

II- 20 PR0P(2,N) Poisson's ratio, v 


NOTES 

(1) MODEL 1 is a linear library model defined by two (2) positive constants 
(E,v) ; i.e., if NPAR(15) .EQ.l , NPAR(17) is set to "2" by default. MODEL 
1 can be used with linear or nonlinear element groups. Since the mat- 
erial constants are independent of history NPAR(18) is "0" by default. 


notes 

( 1 ) 
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VI. 3/D CONTINUUM ELEMENTS (continued) 


MODEL 2 - 
first card 

linear elastic 

orthotropic 


col umn 

variable 

descripti on 

notes 

1-10 

PROP ( 1 ,N) 

E]] Young's modulus in the " 1 " 
principle direction 

( 1 ) 

11-20 

PR0P(2,N) 

E 22 Young's modulus in the "2" 
principle direction 


21-30 

PR0P(3,N) 

E , 3 Young's modulus in the "3" 
principle direction 


31-40 

PRQP(4,N) 

G-j 2 shear modulus 


41-50 

PR0P(5,N) 

G -] 3 shear modulus 


51-60 

PR0P(6,N) 

G 23 shear modulus 


61-70 

PR0P(6,N) 

v-ios Poisson's ratio, contraction in 
11 1 * direction due to stress in " 2 " 
direction 


71-80 

PR0P(7,N) 

v -| 3 Poisson's ratio 


second card 



1-10 

PR0P(8,N) 

\>23 Poisson's ratio 



NOTES 

(1) Model 2 is linear (i.e., NPAR(18).EQ,0 by default). The model 
requires nine (9) property constants, i.e., NPAR (17).EQ.9 by 
default. The constants (E-p, E 22 , . .., V 23 ) are defined in an 
orthogonal principle material coordinate system (11, 22, 33). 
These principle material axes are in turn defined relative to 
an orthorgal coordinate system (a,b,c) associated with each 
element as shown below. 


81 



VI. 3/D CONTINUUM ELEMENTS (continued) 



The "a" direction corresponds to the direction of a vector connecting 
element nodes 6 and 7, "b" is normal to “a" and in> the plane described 
by nodes 5-6-7, and n c" is normal to "a" and "b" such that c = axft. 

The angles a (ALPHA), 3 (BETA) and y (GAMMA) are read individually for 
each element (Section VI. 4) or substructure (Section VII. 2). 

The angles a, 3, and y defining the (1, 2, 3) axes relative to the 
(a, b, c) axes may be defined through a sequence of three axes rotations. 
Assume first that the two coordinate systems coincide. Now rotate the 
(a, b, c) system about the c axis with a degrees to obtain an (a'-b 1 - 
c') system (i.e., a 1 -b 1 is still in the a-b plane and a is the angle 
between a and a 1 )- Next rotate the (a'-b'-c 1 ) system about the b‘ axis 
with 3 degrees to obtain an (a"-b"-c‘ ! ) system (i.e., a"-b'' is now in 
the a'-c' plane and 3 is the angle between c' and c"). Finally rotate 
the (a"-b"~c") system about the a" axis with y degrees to obtain an 
(a"'-b'"-c '") system and, equivalently, the (1-2-3) principal material 
axes. The angles a and y are positive as defined by the right-hand 
rule while 3 is positive using the left-hand rule. 
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VI. 3/D CONTINUUM ELEMENTS (continued) 


MODEL 3 - 
first card 

column 

curve description model 
vari able 

* 

description notes 

1-10 

PROP ( 1 ,N ) 

•i 

volume strain at point "1" (1) 

11-20 

PR0P(2,N) 

2 

Gy f « # • 


21-30 

PR0P(3,N) 

e V 


31-40 

PR0P(4 S N) 

e v 


41-50 

PR0P(5,N) 

e v 


51-60 

PROP (6 »N) 

•5 


61-70 

PR0P(7,N) 

K 1 

k ld 

loading bulk modulus at point"!" 

71-80 

PR0P(8,N) 

K 2 

‘H.D* ' ’ 

• 

second card 




1-10 

PR0P(9,N) 

k ld 


11-20 

PR0P(10,N) 

k ld 


21-30 

P ROP (11 »N) 

k ld 


31-40 

PR0P(12,N) 

k ld 


41-50 

PR0P(13,N) 

r 1 

'UN 

unloading bulk modulus at point "1" 

51-60 

PROP ( 1 4 »N) 

K 2 


61-70 

PR0P(15,N) 

Hin 


71-80 

PROP (16 ,N) 

K 4 

^LIN 


third card 




1-10 

11-20 

PROP ( 1 7, N) 
PROP ( 1 8 iN ) 

K5 

UN 

k un 

(2) 

21-30 

PR0P(19,N) 

g ld 

loading shear modulus at point'"!" 

31-40 

PR0P(20,N) 

r 2 

g ld .. 


41-50 

P ROP (21 ,N) 

r 3 

0 LD 
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VI. 3/D CONTINUUM ELEMENTS (continued) 


column variable description notes 

51-60 PROP (22 ,N) gJ d 

61-70 PR0P(23,N) Gj^ 

71-80 PR0P(24,N) g£ d 


NOTES 

(1) MODEL 3 is a nonlinear library material model for which NPAR(17) .EQ.24, 
be default. To characterize history, thirteen (13) variables— namely, 
the previously calculated stress and strains, and the maximum volumetric 
strain ever reached in loading— need be saved for each integration 
point. Therefore, NPAR(18) is set to "13" by default. 

All six input points must be de fined, as ind icated belowj 



where "J" is the J-th table input point. Linear interpolation is used 
to define the loading and unloading bulk moduli, X K LD> K un , and the 
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VI. 3/D CONTINUUM ELEMENTS (continued) 


loading shear modulus, tg , at time t between input points. If the 

volumetric strain during U time or load step exceeds e^, a diagnostic 
mesage results. 

(2) The unloading shear modulus at time t is calculated as 



^un/WXd. 


NOTATION 


For the notation used in the description of the nonlinear element models 
refer to the report, "Static and Dynamic Geometric and Material Non- 
linear Analysis," SESM Report No. 74-4, Feb. 1974 (Ref. 3). 
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VI. 3/D CONTINUUM ELEMENTS (continued) 

MODEL 4 - Combined isotropic-kinematic hardening, von Mises yield condition. 


CARD 1 


CARDS 

2,3,4 


column 

vari able 

descripti on 

notes 

1-10 

PROP ( 1 ,N) 

Young's modulus, E 

(1) 

11-20 

PR0P(2,N) 

Poisson's ratio, v 

(1) 

21-30 

PR0P(3,N) 

dumber of points in piecewise linear 
stress Strain Curve, RNSSC 

(2) 

31-40 

PR0P(4,N) 

Number of points in piecewise yield 
surface size vs equivalent plastic 
curve, RNSIZE 

(3) 

41-50 

PR0P(5,N) 

Type of hardening rule being used,TPLA$ 

(4) 

51-60 

PR0P(6,N) 

Ratio of isotropic hardening to 
kinematic hardening, RITK 

(5) 

61-70 

PR0P(7,N) 

Value of stress_at first point in stress 
strain curve, o-j 

(6) 

71-80 

PR0P(8,N) 

Value of total strain at first point in 
stress strain curve , 

(6) 

1-10 

PR0P(9,N) 

Value of stress at second point in stress 
strain curve, 


11-20 

PR0P( 10 ,N) 

Value of strain at second point in stress 
strain curve, e 



PROP (5+2* 

RNSSC,N) Value of stress at last point in stress 
strain curve , 

PROP (6+2* 

RNSSC,N) Value of strain at last point in stress 

strain curve, F RNSSC 

PROP (7+2* 

RNSSC,N) Value of yield surface size at first point (7) 
in yield surface size vs equivalent plastic 
strain curve , k-j , 

PROP (8+2* 

RNSSC,N) Value of equivalent plastic strain at first 
point in yield surface size vs equivalent 
plastic strain curve, e] 
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VI. 3/D CONTINUUM ELEMENTS (continued) 


MODEL 4 (continued) 


column variable description notes 


PROP (5+2* 

RNSSC+2*RNSIZE 

,N) Value of yield surface size at last 

point in yield surface size vs 
equivalent plastic strain curve, 

PROP (6+2* 

RNSSC+2*RNSIZE 

,N) Value of equivalent plastic strain at 

last point of yield surface size vs 
equivalent plastic strain curve, ej^SIZE 

NOTES 

(1) MODEL 4 is nonlinear library material model for which NPAR(17) .EQ.26 
by default. NPAR(18) .EQ.21 , since 21 history dependent variables 
must be stored for each integration point. 

The Young's modulus and Poisson's ratio define the initial elastic 
behavior of the material. 

(2) RNSSC is the number of points in a piecewise linear stress strain 
curve, and includes the end point of the last linear segment, but 
not the point STRESS^O, STRAIN=0. Maximum allowable value is 
RNSSC=5, the minimum value is RNSSC=2. 

(3) RNSIZE is the number of points in a piecewise linear yield surface 
size vs equivalent plasti o strain curve, and includes both the 
end point of the last linear segment as well as the point STRESS^ 
yield stress, EQPLAS=0. RNSIZE=0 for TPLAS.NE.3. 

(4) TPLAS EQ.O Isotropic hardening 


EQ.l Kinematic hardening (Prager-Ziegler) 

EQ.2 Combined isotropic, kinematic 

hardening, constant ratio 

EQ.3 Combined isotropic, kinematic 

hardening, variable ratio 
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VI. 3/D CONTINUUM ELEMENTS (continued) 


(5) RITK is the ratio of isotropic to kinematic hardening for TPLAS.EQ.2, 
with RITK. EQ . 1 being pure isotropic, RITK.EQ.O pure kinematic. If 
TPLAS.NE.2, RITK is ignored. 

(6) Values of STRESS and STRAIN are input in pairs for the piecesise 
linear stress strain curve using the format 8F10.0. A maximum 
of five (5) pairs of data is allowed. Input RNSSC pairs of data. 

The input stress-strain curve must be monotonicly increasing such 
that PR0P(5+2*I).LE.PR0P(5+2*(I+1)) and PR0P(6+2*I).LT.PR0P(6+2*(I+1)) 
for any I. See Fig. 9 also. 

(7) Values of YSSIZE (Yield Surface SIZE) and EQPLAS (Equivalent PLAStic 
strain) are input in pairs for the piecewise linear yield surface 
size vs equivalent plastic strain curve. Data for this curve need 
not be input unless TPLAS.EQ.3. For TPLAS.LT.3, a compatable curve 
is generated internally in the program, appropriate for the hardening 
rule indicated by TPLAS. See Fig. 9 also. 

NOTE: Regardless of the value of RNSSC and RNSIZE, four cards will 

be read on 8F10.0 format to generate the PROP array. Blank entries 
for a property value will be read and ignored after the last non-zero 
entry. FOUR CARDS MUST BE SUPPLIED. 
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VI. 3/D CONTINUUM ELEMENTS (continued) 

MODEL 5 - mechanical sublayer model, von Mises yield condition 


column 

vari able 

descri ption 



notes 

1-10 

PROP ( 1 ,N) 

Young's modulus, E 



(1) 

11-20 

PR0P(2,N) 

Poisson's ratio, v 



(1) 

21-30 

PR0P(3,N) 

Number of points in piecewise 
Stress Strain Curve, RNSSC 

linear 

(2) 

31-40 

PR0P(4,N) 

Value of stress at first point 
strain curve, ^ 

in 

stress 

(2) 

41-50 

PR0P(5,N) 

Value of strain at first point 
strain curve , F, 

in 

stress 

(2) 


PR0P(8,N) Value of stress at third point in stress 
strain curve 

\RD 2 PROP (9 S N) Value of strain at third point in stress 

strain curve 


PR0P(2+2*RNSCC 

,N) Value of stress at last point in stress 

strain curve, c Rnssc 

P0RP( 3+2*RNSCC 

,N) Value of Strain: at last point in stress 

strain curve, e RNSSC 


NOTES 

(1) MODEL 5 is a nonlinear library material model for which NPAR(17).EQ.ll 
by default. NPAR(18) -EQ. 30, since 30 history dependent variables 
must be stored for each integration point. 

The Young's modulus and Poisson's ratio define the initial elastic 
behavior of the material. 
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VI. 3/D CONTINUUM ELEMENTS (continued) 

MODEL 5 (continued) 

NOTES 

(2) RNSSC is the number of points in a piecewise linear stress strain 
curve, and includes the end point of the last linear segment (unless 
that segment is perfectly plastic), but not the point $TRESS=0, 
STRAIN s 0. The maximum allowable value is RNSSC=4.0. If RNSSC. LE. 2.0, 
a blank card B must still be supplied. It is also required that the 
input stress strain curve be monotonicly increasing so that PROP (2+2* 
I) .LT. PROP ( 2+2 ( 1+1 ) ) and PR0P(3+2*I).LT.PR0P(3+2*(I+1)). Note that 
the inequality reads .LT., not .LE. See Fig. 10 also. 
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VI. 3/D CONTINUUM ELEMENTS (continued) 

MODEL 6 - Incompressible Nonlinear elastic, Mooney-Ri vlin material 


column 

variable 



description 


notes 

1-10 

PROP ( 1 ,N) 

c i 

C 1 

and C 0 are material 

constants 

0) 


PR0P(2,N) 

in 

the Mooney-Ri vlin expression 


11-20 

C 2 

for strain energy. 



NOTES 







(1) MODEL 

6 is a nonlinear model 

for which NPAR (17) 

is set to 

"a". Since 


the material is hyperelastic, no history-dependent variables need be 
stored at integration points and NPAR(18) is set to "0". C] and C 2 are 
material constants in the expression for strain energy, 

U = C-, ( I-j ~3) + C 2 (I 2 -3) 

where I, and I 2 are the first and second strain invariants, respectively. 
A more complete theoretical description of the model is given in Ref. 5, 
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VI. 3/D CONTINUUM ELEMENTS (continued) 
MODEL 7 - user-supplied material model 


col umn 

vari able 

description 

notes 

1-10 

PROP ( 1 ,N) 


0) 

11-20 

PR0P(2,N) 



21-30 

PR0P(3,N) 




etc. 

NOTES 

(1) MODEL 7 is a non-standard model (i.e., is not in the material model 
library) and requires that property data be read 

( (PR0P(J ,N) ,J=1 ,NPAR(1 7) ,N=1 ,NPAR(16) ) 

where NPAR( 17) is user-supplied (i.e., no default value exists). 
Subroutines which allocate storage locations (in "WA") for element 
history (if the model is nonlinear) and which calculate the incre- 
mental stress-strain law and -current stresses must be supplied. 
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VI. 3/D CONTINUUM ELEMENTS (continued) 

VI. 3. Stress Output Table Cards (1615)* 

Skip this section if stresses are to be printed at integration points, 
i.e., if NPAR( T 3) .EQ.O; otherwise supply NPAR(13) cards. 


column 

vari able 



description 

notes 

1-5 

ITABLE (N,l) 

Stress 

output 

location 

1 

(1) 

6-10 

ITABLE(N,2) 

Stress 

output 

locati on 

2 


11-15 

ITABLE(N,3) 

Stress 

output 

location 

3 


16-20 

ITABLE(N ,4) 

Stress 

output 

location 

4 


21-25 

ITABLE (N ,5) 

Stress 

output 

location 

5 


26-30 

ITABLE(N ,6) 

Stress 

output 

location 

6 


31-35 

ITABLE (N, 7) 

Stress 

output 

location 

7 


36-40 , 

ITABLE (N, 8) 

Stress 

output 

location 

8 


41-45 

ITABLE (N, 9) 

Stress 

output 

location 

9 


46-50 

ITABLE(N,10) 

Stress 

output 

location 

10 


51-55 

ITABLE (N ,11) 

Stress 

output 

location 

11 


56-60 

ITABLE(N,12) 

Stress 

output 

location 

12 


61-65 

ITABLE(N,13) 

Stress ■ 

'output 

location 

13 


66-70 

ITABLE(N,14) 

Stress output 

location 

14 


/ 1-75 

ITABLE(N ,15) 

Stress 

output 

location 

15 


76-80 

ITABLE(N ,16) 

Stress 

output 

location 

16 



NOTES 

(!) Stress tables are defined to provide flexibility in element stress 

output requests. Each element can refer to one of the tables defined, 
and the element stresses are then calculated at the points specified 
in that table. Figure 15 defines the locations of the 27 possible 
points within an element where stresses may be printed. 

Any one table may contain a maximum of sixteen (16) stress output 
points. The first "0" entry in a table will terminate that table. For 
example, if ITABLE (N, 1 ),. EQ.21 and ITABLE(N,2).EQ.O, then stresses will 
be printed only at point 21 (the centroid of the element) whenever this 
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VI. 3/D CONTINUUM ELEMENTS (continued) 


NOTES 


stress table is referred to. 

Stress output tables are used only with linear material models (NPAR(15) 
.EQ.l or l). For elements using nonlinear material models (NPAR(15) . GE. 3) , 
stresses can be printed at integration points only. 
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VI. 3/D CONTINUUM ELEMENTS (continued) 


VI. 4. Element Data Cards* 

Three cards must be prepared for each element that appears in the 
input: 

CARD 1 (715, 3F10.0) 


column variable 

I- 5 M 

6-10 IELD 

II- 15 IELX 

16-20 IPS 

21-25 MTYP 

26-30 1ST 

31-35 KG 

36-45 ALPHA 

46-55 BETA 

56-65 GAMMA 


description notes 

3/D CONTINUUM element number (1) 

GE.l and LE.NPAR(2) 

Number of nodes to be used in des- (2) 

cribing the element's displacement 

field 

EQ.O - default set to NPAR(7) 

LE.NPAR(7) 

Number of nodes to be used in the (3) 

description of the element's geometry 

EQ.O - default set to IELD 

GE.8 and LE.IELD 

EQ.IELD ■* isoparametric element 

LT.IELD -* subparametric element 

Number of the stress output location (4) 

table to be used for this element (if 
a linear material model is used) or flag 
for printing stresses at element i nte- 
gration points 

EQ.O - no stress output for this element 

Material property set number assigned to 

this element 

GE.l and LE. NPAR(16) 

Flag indicating that the stiffness and (5) 
mass matrices for this element are the 
same as those for the preceding element 
EQ.O - no 
EQ.l - yes 

Node number increment for element (6) 

data generation 

EQ.O - default set to"l" 

Angle a (degrees) used to define princi- (7) 
pal material axes for Model 2 (orthotropic) 

Angle B (degree) used to define principal 
material axes for Model 2 (orthotropic) 

Angle y (degrees) used to define principal 
material axes for Model 2 (orthotropic) 
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VI. 3/D CONTINUUM ELEMENTS (continued) 

CARD 2 (815) 

column variable description notes 


1-5 

N0D(1) 

Global node 
point 1 

6-10 

NOD ( 2 ) 

Global node 
point 2 

11-15 

N0D(3) 

Global node 
point 3 

16-20 

N0D(4) 

Global node 
point 4 

21-25 

N0D(5) 

Global node 
point 5 

26-30 

N0D(6) 

Global node 
point 6 

31-35 

N0D(7) 

Global node 
point 7 

36-40 

N0D(8) 

Global node 
point 8 


CARD 3 (1315) 


column 

variable 


1-5 

N0D(9) 

Global node 
point 9 

6-10 

■NOD(IO) 

Global node 
point 10 

11-15 

NOD(ll) 

Global node 
point 11 

16-20 

N0D(12) 

Global node 
point 12 

21-25 

N0D( 1 3) 

Global node 
point 13 

26-30 

N0D(14) 

Global node 
point 14 

31-35 

N0D(15) 

Global node 
point 15 


number of element nodal (2) 

number of element nodal 

number of element nodal 

number of element nodal 

number of element nodal 

number of element nodal 

number of element nodal 

number of element nodal 

description notes 

number of element nodal (2) 

number of element nodal 

number of element nodal 

number of element nodal 

number of element nodal 

number of element nodal 

number of element nodal 
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VI. 3/D CONTINUUM ELEMENTS (continued) 


column 

variable 


36-40 

N0D(16) 

Global node 
point 16 

41-45 

N0D( 1 7} 

Global node 
point 17 

46-50 

N0D( 18) 

Global node 
point 18 

51-55 

N0D( 19) 

Global node 
point 19 

56-60 

N0D(20) 

Global node 
point 20 

61-65 

NOD (21 ) 

Global node 
point 21 

NOTES 




description notes 

number of element nodal 

number of element nodal 

number of element nodal 

number of element nodal 

number of element nodal 

number of element nodal 


(1) Element cards must be input in ascending element number order 
beginning with one (1) and ending with NPAR(2). Repetition of 
element numbers is illegal, but element cards may be omitted, 
and missing element data are generated according to the procedure 
described in note (6). 


(2) IELD is a count of the node numbers actually posted on Cards 2 
and 3 which must immediately follow Card 1. IELD must be_at 
least eight (8), but must be less than or equal to the limit 
NPAR(7) which was given on the element group control card. Section 
VI. 1. Element displacements are assigned at the IELD non-zero 
nodes, and thus the order of the element matrices is three (i.e., 
translations X,Y,Z) times IELD. The eight corner nodes of the 
hexahedron must be input, but nodes 9 to 21 are optional, and any 
or all of these optional nodes may be used to describe the element's 
displacement field. However, all 21 entries for N0D(i ) are read 
from element data cards 2 and 3; if IELD. LT. 21 the particular node 
locations not used in this element must be input as zero (0) in 
NOD(i)- Figure 16 defines the input sequence that must be observed 
for element input. 
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VI. 3/D CONTINUUM ELEMENTS (continued) 


NOTES 


For example, the 10-node element (IELD.EQ.10) shown below 



is defined by 


N0D(1) = [X X X X X X H 0 0 X 0 0 0 X 0 0 0 0 0 O], 

where the nonzero entries (X) are the global mesh node numbers 

(•Section III) of the 10 nodes. 

(3) When element edges are straight it is unnecessary computationally 

to include side nodes in the numerical evaluation of coordiante 
derivatives, the Jacobian matrix, etc., and since regular element 
shapes are common, an option has been included to use fewer nodes 
in these geometric calculations than are used to describe element 
displacements. The first IELX nonzero nodes posted on Cards 2 and 3 
are used to evaluate those parameters which pertain to element 

geometry only. IELX must be at least eight (8), and if omitted is 

re-set to IELD. A common application might be a 20 node element 
(i.e., IELD. EQ. 20) with straight edges in which case IELX would be 
entered as "8". 

(4) If IPS.EQ.O no stress output will be provided for this element. A 
nonzero entry for the element stress printout flag "IPS" has a dif- 
ferent meaning when used with nonlinear material models than it does 
when used with linear material models: 

a) For a linear material model (MODEL . EQ . 1 or 2) "IPS" specifies 
the stress output table (Section VI. 3) to be used for this 
element. If no stress output tables were input (i.e., NPAR(13). 
EQ.O) then IPS - EQ. 1 specifies that stresses are to be printed 
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VI. 3/D CONTINUUM ELEMENTS (continued) 


NOTES 

at all integration points for this element. 

b) For a nonlinear material model (MODEL. GE. 3) IPS.EQ.l specifies 
that stresses are to be printed at all integration points for 
this element. 

(5) The flag 1ST allows the program to bypass stiffness and mass matrix 
calculations providing the current element is identical to the 
preceding element; i.e., the preceding and current elements are id- 
entical except for a rigid body translation. If IST.EQ.O, new 
matrices are computed for the current element. 

(6) When element cards are omitted, element data are generated automatically 
as follows, 

a) All data on Card 1 for generated elements is taken to be the 
same as that given on the first element card in the sequence 

b) Nonzero node numbers (given on Cards 2 and 3 for the first element) 
are incremented by the value "KG" (on the first element's Card 1) 

as element generation progresses; zero (or blank) node number entries 
are generated as zeroes. 


The last element cannot be generated and must be included. 


(7) The angles a, $ and y are required only for Model 2. They are defined 
in the notes for Model 2. 
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VII. NODE AND ELEMENT INPUT BY SUBSTURCTURES 


If the nodal coordinates, material property cards and element-node 
connectivity are to be input by using the substructure option (MG0=1), 
then this section (VII)is input in place of Section III-VI. 

In brief, a complete structure is described as a collection of sub- 
structures. Each substructure is described by specifying its corner nodes 
and up to two midside nodes for each side, the number of mesh points to be 
generated in each coordinate direction, and various gradients to control the 
mesh spacing. Each substructure may be thought of as a single isoparametric 
element which is to be broken up into a number of smaller elements. The 
program then uses isoparametric shape functions to interpolate nodal co- 
ordinates and to generate element-node connectivity for all elements within 
the substructure. It is assumed that each element within the substructure 
has identical nonlinear material parameters, the same number of nodes, etc. 

As an example, consider the pressure vessel shown below: 
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The above structure may be conviently described by five substructures as 
shown. The required input would consist of the corner and midside nodes of 
each substructure plus information describing the number of mesh points in 
each of the coordinate directions (local element coordinates, say r and s) 
and the gradients controlling mesh size. 

Enter one set of the cards in Section VII for each substructure. Linear 
substructures (small strain, elastic) must be input first followed by 
nonlinear substructures. Substructures must be input in sequence so that 
they attach to the right, top, or back (3-D) of a previously input sub- 
structure. The side of a substructure which attaches to a previously input 

substructure must have identical corner nodal coordinates and number of 

/ 

nodes along the matching boundary, i.e., the entire side must be matched. 
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VII. 1 Control Card (1015, 6F5.2) 


column variable description notes 

I- 5 NPAR(l) Element Type (1) 

EQ.l - truss 

EQ.2 - 2-D continuum 

EQ.3 - 3-D continuum 

6-10 NPAR(3) Geometric nonlinear analysis type (2) 

EQ.O - small strain (linear) 

EQ.l - large strian (total Lagrangian) 

EQ.2 - large strain (updated Lagrangian) 

II- 15 NPAR(5) Element type code (2-D only) 

EQ.O - axi symmetric 

EQ.l - plane strain 

EQ.2 - plane stress 

16-20 NPAR(IO) Integration order (for r-s plane) (3) 

EQ.O - default reset to 2 
GE.2 and LE.4 

21-25 NPAR(ll) Integration order for t direction (4) 

(3-D only) 

EQ.O - default reset to 2 
GE.2 and LE . 4 

26-30 NPAR(15) Material model number (see table) (5) 

31-35 ITYPE Number of boundary nodes used to (6) 

describe this substructure: 

For 2-D; ITYPE=4 (linear), 8 (quadratic) 
or 12 (cubic) 

For 3-D, ITYPE=8 (linear, 20 (quadratic) 
or 32 (cubic) 


36-40 

I MESH 

Number of mesh points to be generated in X (7) 
direction (3-D only) 

41-45 

JMESH 

Number of mesh points to be generated in Y 
di recti on 

46-50 

KMESH 

Number of mesh points to be generated in Z 
direction 

51-55 

G1 

Gradient used to determine mesh spacing (8) 

from node 1 to node 2 of substructure 
EQ.0.0 - default reset to 1.0 

56-60 

G2 

Gradient used to determine mesh spacing 
along top, front edge from node 1 to node 4 
EQ.0.0 - default reset to 1.0 
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VII. 1 CONTROL CARD (1015, 6F5.2) (continued) 


column 

variable 

description 

notes 

61-65 

G3 

Gradient used to determine mesh spacing 
along top, left edge from node 4 to node 
EQ.0.0 - default reset to 1.0 

3 

66-70 

G4 

Gradient used to determine mesh spacing 
along top, rear edge from node 2 to node 
EQ.0.0 - default reset to 1.0 

3 

71-75 

G5 

Gradient used to determine mesh spacing 
along edge from node 1 to node 5 
EQ.0.0 - default reset to 1.0 


76-80 

G6 

Gradient used to determine mesh spacing 



along edge from node 4 to node 8 
EQ.0.0 - default reset to 1.0 


NOTES 

(1) All elements that are generated within this substructure will be of 
the same type, i.e., will have the same parameters for NPAR(l), 

NPAR( 3) , NPAR(5), NPAR(IO), NPAR(ll), NPAR(15), ITYPE , the same 
material model properties, etc. 

(2) NPAR(3) determines whether geometric nonlinearities (large dis- 
placements and large strains) are to be included in the analysis. 

(3) For rectangular elements, a Guassian numerical order of 2 is generally 
sufficient. For distorted elements (curved boundary or skewed sides), 
a higher integration order is needed to accurately represent the 
volume integration. It should be noted that the computation time 

for stiffness and force vectors is strongly dependent on the inte- 
gration order. Thus, 3 by 3 integration requires more than twice 
the computer time that 2 by 2 integration requires. 

(4) This entry is required only if this substructure represents 3-D 
elements (NPAR(l) .EQ. 3) . For some applications such as plates or 
shells where the strain varies approximately linearly through the 
"thickness", the integration order through the thickness may be 
less than over the surface (r-s plane) of the shell. 

(5) The available material models are tabulated in Table 2. 

(6) The substructure can be thought of as a large isoparametric element 
whose substructure geometry is described by boundary nodes in exactly 
the same manner as a single element is described. The elements gen- 
erated within a substructure are given the same number of nodes as the 
parent substructure. Hence, a linear substructure (4 nodes for the 
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VII. 1 CONTROL CARD (1015, 6F5.2) (continued) 


NOTES 


the 2-D case and 8 nodes for the 3-D case) should be used only if it 
is desired to generate linear elements within the substructure. 

(7) IMESH, JMESH, and KMESH specify the number of nodal points to be 
generated in the r,s and t directions, respectively. If EQ.O, no 
mesh is generated and the substructure is treated as a single element. 
When a substructure boundary is to match or mate together with a 
previously input substructure, the number of generated nodes must be 
identical on the adjacent boundaries. Subsequently input substructures 
must mate on the top, right or back boundaries of previously input 
substructures. See Figs. 17 and 18. 


(8) The gradient parameters SI through G6 are used to determine size 
variation of the generated mesh (node points). If the gradient is 
greater than unity, the element size will decrease away from the 
origin. If the gradiant is less than unity, the element size will 
increase away from the origin. Examples of the generated mesh for 
various gradients are shown in Figure 18. 
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BOTTOM : 3-4 
RIGHT: 4-1 


Fig. 17 ‘ Gradient Parameters for Two-Dimensional Substructure 
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Fig. 19 Example of Generated Two-Dimensional Mesh 
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VII. 2 Material Property Control Card (5F10.0, 4I5)-0MIT THIS CARD FOR A 
TRUSS ELEMENT GROUP (NPAR(1 )=1 ) . 


column variable 


description 


notes 


I- 10 THICK 

II- 20 DENSE 

21-30 ALPHA or 

BETA 


31-40 BETA or 
STIFAC 


41-50 GAMMA or 
SHE FAC 


51-55 IPS 


56-60 NPAR(16) 

61-65 NPAR(1 7) 


Thickness of element (plane stress only) (1) 

Mass density (required only in dynamic (2) 

analyses or if body forces are included 


Dual meaning entry depending upon material (3) 
model : 

3-D, Model 2 - ALPHA - angle a (degrees) 

used to define 
principal material 
axes 

2- D, Model 2 - BETA - orthotropic mat- (4) 

erial angle 8 

Dual meaning entry depending upon material 
model : 

3- D, Model 2 - BETA - angle 8 (degrees) (3) 

used to define prin- 
cipal material axes 

2- D, Model 4,5 - STIFAC - Stiffness re- (5) 

duction factor 
after fracture 

Dual meaning entry depending upon material 
model : 

3- D, Model 2 - GAMMA - Angle y (degrees) (3) 

used to define 
principal material 
axes 

2-D, Model 4,5 - SHEFAC - Shear .reducti on (5) 

factor after 
fracture 

Stress print out flag. (6) 

For linear material models, refers to 
number of stress output table to be used 
for all elements in this group. For non- 
linear material models, flag for printing 
stresses at integration points. 

EQ.O - no stress output for this substructure 

Number of material property sets (must (7) 
be 1 for MG0=1 option) 

Number of constants per property set (8) 

(required only for used supplied material 
models) 
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VII. 2 Material Property Control Card (5F10.0, 415) (continuted) 


col umn 

variable 

description 

notes 

66-70 

NPAR( 1 8) 

Dimension of storage array required 
for element history-dependent 
variables (required only for user 
supplied material models). 

(8) 


NOTES 

(1) The element thickness is required for plane stress elements. Plane 
strain elements are assumed to have a thickness of one unit and 

axi symmetric elements have a thickness of one radian. 

(2) The mass density is used in dynamic analyses to compute the mass 
matrix. It is also used in 2-D analyses when body forces are in- 
cluded. The mass density must be input In proper units; no con- 
version factors are applied by the program. 

(3) The angles a, B and y are required only for 3-D, Model 2 only. 

They are defined in the notes for this material model (Section 
VI. 2. b. Model 2). 

(4) The orthotropic material angle B is required for 2-D, Model 2 only. 

It is defined further in the notes for this material model (Section 
V.2.b, Model 2). 

(5) STIFAC and SHEFAC are assumed material parameters used in conjunction 
with 2-D analyses with M0DEL=4 or 5. See the notes for these parameters 
in Section V.2.b (MODEL 4). 

(6) If IPS.Eq.O, no stress output will be provided for the elements in this 
substructure. A nonzero entry has a different meaning when used with 
linear or nonlinear material models as follows: 

For a linear material model (MODEL. EQ.l or 2), IPS specifies the stress 
output table to be used for elements in this substructure. If no stress 
output tables were input (i.e., NPAR(13).EQ.O) , then IPS. EQ.l means that 
stressea are to be printed at all integration points of each element. 

For a nonlinear material model (MODEL. GE. 3) , IPS, EQ.l specifies that 
stressea are to be printed at all integration points of each element. 

(7) It is assumed that each element in the substructure has identical mat- 
erial properties and therefore NPAR(16) = 1. 

(8) For standard material models, in the program, these entries may be 
omitted. For user supplied models, the notes regarding these entries in 
Section V.l and VI. 1 should be consulted. 
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VII. 3 Material Property Cards (and Stress Output Tables) 


Material property .cards are required for the particular material 
model (NPAR( 15) ) being^uti lized in this substructure. The required 
input is identical to that described in Sections IV. 2, V.2, and VI. 2. 
Depending upon the type of element, the material property card and 
stress output table input described in the following sections should 
be input: 


1- D: Sec. IV. 2. b. Models 1-2 

2- D: Sec. V.2.b and V.3, Models 1-10 

3- D: Sec. VI. 2. b and VI. 3, Models 1-6 
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VII. 4. Nodal Coordinates of Substructure (II, 13, 311, 3F20.0) 


The following card is repeated for as many nodes as are being used 
to describe this substructure (ITYPE on Card VII.l). 


column 

variable 

descri ption 

notes 

1 

IT 

Coordinate system used to describe 
this node 

EQ.O - cartesian (X,Y,Z) 

EQ.l - X-cylindrical 

(1) 

2-4 

N 

Node number 

(2) 

5 

I XT RAN 

X-translati on boundary condition code 
EQ.O - free 
EQ.l - fixed 

(3) 

6 

IYTRAN 

Y-translation boundary condition code 


7 

IZTRAN 

Z-translation boundary condition code 


11-30 

X (N ) 

X coordinate; or Z coordinate if IT. EQ.l 
(used only in 3-D) 

(4) 

31-50 

Y(N) 

Y coordinate; or R coordinate if IT. EQ.l 


51-70 

Z(N) 

Z coordinate; or e coordinate if IT. EQ.l 



NOTES 

(1) Nodal coordinates may be input in either the global, cartesian or 
X-cylindrical coordinate systems. A complete description of the 
X-cylindri cal coordinate system is given in Note 1 following 
Section III.l. 

(2) Node numbers must range from 1 to ITYPE (Sec. VII.l) but need not be 
input in ascending order. Figures 17 and 18 define the node numbering 
conventions that must be followed. 

(3) Boundary condition codes are used to delete (fix) certain degrees- 
of-freedom. If the same boundary conditions are specified for all 
the nodes on one edge of an element, then these boundary conditions 
are applied to all nodes generated by the program on that edge. If, 
for example, boundary conditions are applied only to one node on 

a particular edge, then these boundary conditions are retained for 
that node but all other generated nodes on that edge are given free 
boundary conditions. 
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VII. 4. Nodal Coordinates of Substructure (continued) 


NOTES 

(4) The nodal coordinates should be input in the proper coordinate 
system as prescribed by the coordinate system flag IT in CC 1. 
The origin of the global, cartesian and cylindrical coordinate 
systems are the same. It should be noted that the boundary 
condition codes always refer to degrees -of- freedom in the global 
X-Y-Z system even if the nodal coordinates were given in the 
cylindrical coordinate system. 
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VIII. LOAD CURVES (continued) 


Hence, the definition of the applied pressure (a constant) and the 
magnitude of the load curves points is arbitrary. The only require- 
ment is that the product of the constant pressure and the load curve 
multiplier evaluated at a particular time yield the correct magnitude 
that is to be applied at that time. 
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IX. NONCONSERVATIVE PRESSURE LOADS 

Omit these cards if NPRESN.EQ.O; otherwise, input NPRESN cards. 

Nonconservative pressure loads are assumed to be functions of the de- 
formation, i.e., normal pressure loads remain normal to the deformed 
surface, etc. Nonconservative pressure loads can be applied only to 
elements for which large deformations (strains) are being considered. 

The same nonconservative pressure load, i.e., one card, may be used to. 
describe loads being applied to one element, to a set of elements within 
a group or to an entire group (substructure). 


col umn 

variable 

descripti on 

notes 

1-5 

IGROUP 

Substructure or element group 
to which load is applied 

(1) 

6-10 

IELM1 

First element to which load is 
applied (ignored if MG0=1) 

(2) 

11-25 

IELM2 

Last element to which load is 
applied (ignored if MG0=1) 

(2) 

18-25 

BONDRY 

Boundary of substructure of element 
to which load is applied. Entry 
must be right justified; 

(3) 


TOP 

LEFT 

BOTTOM 

RIGHT 

™^ T |only in 3-D elements 


26-30 

ICUR 

Load curve number according to 
which the load is applied 

(4) 

31-40 

PRESN 

Normal pressure (positive if 
tensile) 

(5) 

41-50 

PRESS 

Shear pressure (positive according 
to Timoshenko notation) - shear 
pressure not allowed in 3-D elements 

(5) 


NOTES 

(1) IGROUP must be a valid group number, i.e., 1 <_ IGROUP <NUMEG. 

The same group number may be repeated if necessary 

(2) IELM1 and IELM2 define an inclusive set of elements to which the 
pressure is to be applied. When describing the mesh by sub- 
structures, the program assumed that all elements in substructures 
IGROUP are loaded, i.e., these parameters are ignored if MG0=1. 
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VIII. LOAD CURVES* 


If no load curves are to be input, omit the cards in Sec. VII. Otherwise, 
input one card VIII. 1 followed by one or more cards VIII. 2; then repeat 
this sequence until data for NLCUR load curves has been input. 


VIII. 1. 

Control card (15) 


column 

vari able 

•description 

notes 

1-5 

NPTS 

Number of points used to input this 
load curve 
GE.2 and LE. 100 


VIII. 2. 

Load Curve 

Multiplier Data (8F10.0) 


column 

variable 

description 

notes 

1-10 

TIMV(l) 

Load (time) at point 1 

(1) 

11-20 

RV(1) 

Load multiplier at point 1 


21-30 

TIMV(2) 

Load (time) at point 2 


31-40 

RV(2) 

Load multiplier at point 2 



71-80 RV{4) 

Next card (if required) 

I- 10 TIMV (5) 

II- 20 RV(5) 


NOTES 

(1) The data pairs [TIMV(i), RV(i)] define load multiplier functions 
which describe the variation of applied loads with time (or.a load 
parameter in static problems). The magnitude of the load at any 
time (load parameter) is obtained by linearly interpolating the function 
and multiplying this times a constant input pressure assigned to an 
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VIII. LOAD CURVES (continued) 
NOTES 


element or a force assigned to a node (as input on cards IX, X, or 
XI). The first time, i.e., TIMV(l) must be equal to zero. The last 
time value for any curve, i.e., TIMV(NPTS) must be greater than or 
equal to the time at the end of the solution, i.e., TIMV(NPTS) >_ 
TSTART + NSTE * DT. 

A typical load curve is shown below. The number 



of points to be input is four. Suppose a certain element has a 
pressure of 500 psi which is applied according to the above load curve. 
At a time (load) of 1.0, the magnitude of the pressure will be 
0.5 x 500 = 250 psi. At a time' (load) of 2.0, the magnitude will be 
0.75 x 500 = 375 psi, etc. Note that the same magnitudes would be 
obtained if the pressure specified for the element were 100 psi and 
the load curve was input as 




IX. NONCONSERVATIVE PRESSURE LOADS (continued) 


NOTES 

(3) The pressure load on this card is applied to the face of an element 
(elements or substructure) as shown inFigs.20 & 21. Note that the 
faces are defined with respect to the element axes (r,s,t) which 
are in turn defined by the element node numbers 

(4) The pressure on this card is applied to the element face with a time 
variation described by load curve ICUR. Note that 1 £ ICUR £ NLCUR. 

(5) The normal and shear pressure loads are positive as shown in Figs. 

20 and 21. 
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BOTTOM : 3 -4 
RIGHT: 4-1 


Fig. 20 Nomenclature for Inputting Pressure Distribution on 
Two-Dimensional Element 
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FRONT 



TOP • i-2-3-4 

LEFT *- 3-4-8 -7 

BOTTOM : 5- 6 -7- 8 
RIGHT: I-5-6-2 

FRONT : 1-4-8-5 


BACK : 2-3-7-S 


Fig. 21 Nomenclature for Inputting Pressure Distribution on 
Three-Dimensional Element 


X. CONSERVATIVE PRESSURE LOADS 

Omit these cards if NPRESC.EQ.O; otherwise, input NPRESC cards. 

Conservative pressure loads are assumed to be independent of deformation. 

The same conservative pressure load, i.e., one card, may be used to 
describe loads being applied to one element, to a set of elements within 
a group, or to an entire group (substructure). 


col umn 

variable 

description 

notes 

1-5 

JGROUP 

Substructure or element group to which 
load is applied 

0) 

6-10 

I ELM! 

First element wo which load is 
applied (ignored if MG0=1) 

(2) 

11-15 

IELM2 

Last element to which load is 
applied (ignored if MG0=1) 

(2) 

18-25 

BONDRY 

Boundary of substructure of element 
to which load is applied. Entry 
must be right justified: 

TOP 

LEFT 

BOTTOM 

RIGHT 

BAc!<T } on1y in 3-0 e1ements 

(3) 

26-30 

ICUR 

Load curve number according to which the 
load is applied 

(4) 

31-40 

PRESN 

Normal pressure (positive if tensile) 

(5) 

*1-50 

PRESS 

Shear pressure (positive according to 
Timoshenko notation) - shear pressure 
not allowed in 3-D elements 

(5) 

NOTES 





(1) JGROUP must be a valid group number, i.e., 1 £ JGROUP NUMEG. 

The same group number may be repeated if necessary. 

(2) IELM1 and IELM2 define an inclusive set of elements to which the 
pressure is to be applied. When describing the mesh by sub- 
structures, the program assumed that all elements in substructures 
IGROUP are loaded, i.e., these parameters are ignored if MG0=1. 

(3) The pressure load on this card is applied to the face of an element 
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X. CONSERVATIVE PRESSURE LOADS (continued) 
NOTES 


(elements or substructure) as shown in Fig. IX. 1. Note that the 
faces are defined with respect to the element axes (r,s,t) which 
are in turn defined by the element node numbers as shown in Fig. IX. 1. 

(4) The pressure on this card is applied to the element face with a time 
variation described by load curve I CUR. Note that 1 < ICUR 5 NCLUR. 

(5) The normal and shear pressure loads are positive as shown in Figs. 20 
and 21 . 
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XI. CONCENTRATED NODAL FORCES 

Omit these cards if NCONFC.EQ.O; otherwise, input NCONFC cards. 

The concentrated nodal forces are assumed to be conservative and 
independent of deformation. 


col umn 

variable 

description 

notes 

1-5 

JGROUP 

Substructure to which load is 
applied (ignored if MG0=0) 

(1) 

6-10 

NOD 

Node to which load is applied. 
For MG0=0, NOD corresponds to 
the global node numbers as input 
in Section III. 

For MG0=1 , NOD corresponds to 
the local node number of the 
substurcture. 

(2) 

11-15 

IDIR 

Direction of load 
EQ.l - X direction 
EQ.2 - Y direction 
EQ.3 - Z direction 

(3) 


16-20 

JLCUR 

Load curve according to which load 
is applied 

(3) 

21-30 

FORCE 

Nodal force 

(4) 


NOTES 

(1) JGROUP must lie in the range 1 £ JGROUP £ NUMEG. The same group 
number may be repeated as necessary. 

( 2 ) The location and direction of the concentrated nodal force is given 

by NOD and IDIR, respectively. For the node input option, MG0=0, NOD is 
a global node number as input in section and 1 £ NOD £ NUMNP. For 
the mesh generation option using substructures, MG0=1, NOD corresponds 
to a local element node number as input in Section VII; valid range 
is 1 £ NOD £ 8 for 2-D elements and l£ NOD £ 20 for 3-D elements. 

The direction of the load is specified by IDIR and is always with 
respect to the global axes direction. Note that for 2-D problems, 
the geometry is defined with respect to the (Y-Z) plane and loads 
may only be applied in the Y and Z direction. 

(3) The concentrated force is to vary with time (or load parameter) 
according to load curve JLCUR. Note that 1 £ OLCUR £ NLCUR. 

(4) The nodal, force is positive in the global (X,Y,Z) directions. 
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XII. RAYLEIGH DAMPING (2F10.0.f 

Omit this card if no Rayleigh damping is being included, IDAMP.EQ.O. 


column 

variable 

description 

notes 

1-10 

ADAMP 

Rayleigh damping coefficient a 

0) 

11-20 

BDAMP 

Rayleigh damping coefficient g 

0) 


NOTES 

(1) Rayleigh damping is defined as C = aM + $K, where a and $ are input 
as above. 

It should be noted that a and g are applied to the linear mass and 
stiffness matrices. 


+ Not allowed in current program version. 
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XIII. CONCENTRATED NODAL MASSES (110, 6F10.0) f 

Omit this section if IMASSN.EQ.O; otherwise, input IMASSN cards as follows. 


column 

variable 

description 

notes 

1-10 

N 

Node number 

(1) 

11-20 

XMASS(l) 

X-di recti on mass 

(2) 

21-30 

XMASS(2) 

Y-di recti on mass 


31-40 

XMASS(3) 

Z- direction mass 


41-50 

XMASS(4) 

X-rotational mass 

(3) 

51-60 

XMASS(5) 

Y-rotational mass 


61-70 

XMASS ( 6 ) 

Z- rotational mass 


NOTES 




(1) Input IMASSN cards; 
accumulates mass at 

node order is not important. Repeating nodes 
that node. The node number corresponds to the 


global node number; hence 1 _< N £ NUMNP. If the substructure option 
is used to generate node numbers, a data check run must be made 
first to determine the global node number assigned by the mesh gen- 
erator. Once the node numbers are known, the node numbers required 
on this card may be input. 


(2) Mass components input for deleted (or non-existent) degrees of 
freedom are ignored by the program without a diagnostic message. 

(3) Rotational degrees of freedom are currently not used. 


+ Not allowed in current program version. 
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XIV. CONCENTRATED NODAL DAMPERS (110, 6F10.0) + 

Omit this section if IDAMPN.EQ.O; otherwise, input IDAMPN cards as follows. 


col umn 

variable 

description 

notes 

1-10 

N 

Node number; 

6E.1 and LE.NUMNP 

(1) 

11-20 

XDAMP (1) 

X-di recti on damper 

(2) 

21-30 

XDAMP(2) 

Y- direction damper 


31-40 

XDAMP ( 3) 

Z-di recti on damper 


41-50 

XDAMP(4) 

X- rotation damper 

(3) 

51-60 

XDAMP (5) 

Y- rotation damper 


61-70 

XDAMP (6) 

Z-rotation damper 


NOTES 

(1) Input 

IDAMPN cards; 

node order is not important. 

Repeating nodes 


accumulates dampers at the node. 

If the substructure option is used to generate node numbers, a data 
check run must be made first to determine the global node numbers 
assigned by the mesh generator. Once the node numbers are known, 
the node numbers required on this card may be input. 

(2) Damper components input for deleted (or non-existent) degrees of 
freedom are ignored by the program without a diagnostic message. 

(3) Rotational degrees of freedom are currently not used. 


•j* 

Not allowed in current program version. 
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XV. INITIAL CONDITIONS 


Initial conditions for all active degrees -of- freedom are defined in this 
section. If this is a restart run (M0DEX.EQ.2) , the cards in this section 
must be input but are ignored by the program. The variable ICON controls 
whether zero initial conditions are automatically generated or non-zero 
initial conditions are read from data. 


XV. 1. Control Flag (15) 


column variable description notes 

1-5 ICON Flag indicating the type of 

conditions to be input initial 
EQ.O - zero initial conditions are 
generated automatically 
EQ.l - non-zero initial cnoditions 
are read from data cards. 


XV. 2. Card Input of Initial Conditions (6E12.6) 

If ICON. EQ.l, the program performs the following read operations: 


READ(5,1000) (DIS(K, K = 1 ,NEQ) 

READ(5,1000) (VEL(K) ,K = 1 ,NEQ) 

READ(5,1000) (ACC(K) ,K = 1 ,NEQ) 

1000 FORMAT (6E12.6) 

where DIS, VEL, and ACC are the global displacement, velocity, and ac- 
celeration vectors, respectively. The variable NEQ is the total number 
of degrees-of-freedom, i.e., six times the total number of nodes minus 
all deleted degrees-of-freedom provided by fixed boundary conditions. 

The list of equation numbers assigned to each degree-of- freedom may be 
obtained from the print-out of a previous run. It may also be determined 
be careful inspection of the node numbers and boundary conditions input 
by the program user. 

For a static solution (IMASS.EQ.O) , the VEL and ACC initial vectors are 
not read from card input. 
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XVI. FREQUENCY SOLUTION DATA CARD (51, FIO.O, I5) + 

Input this card only if a frequency and mode shape solution is to be 
performed, i.e., if IEIG.EQ.l 


col umn 

variable 

description 

notes 

1-5 

NFREQ 

Number of frequencies to be calculated 
EQ.O - default set to 1 

(1) 

6-15 

COFQ 

Cut-off circular frequency (units of 
radians per unit time) 

EQ.O - default set to 1.0E + 08 

(1) 

16-20 

IFPR 

Flag to suppress intermediate print- 
out of solution variables, interim 
informative messages, etc.; 

EQ.O - print-out suppressed 
EQ.l - print-out produced 

(2) 


NOTES 

(1) The lowest "NFREQ" vibration mode shapes and frequencies are evaluated 
evaluated for the structure linearized at time "TSTART". Frequencies 
are extracted in numerically ascending sequence, and unless the "N-th" 
frequency is > COFQ, the program will continue the eigensolution 
until all "NFREQ" modes have been determined. If the N'th frequency 
is > COFQ, where 1 ^N <_ NFREQ, the eigensolution is discontunued, 

and only the lowest "N 1 frequencies and mode shpaes are printed. 

At completion of solution, "physical" error bounds for each eigenvalue 
and eigenvector are calculates using 1 1 (K-a> 2 )<j>| |/ j | K<j>[ | , where 
(a),<|>) is a calculated eigenpair. 

(2) During the course of eigenvalue extraction various solution variables, 
informative diagnostics, etc. are printed if IFPR.EQ.l. This print- 
out is useful in tracing solution behavior. 


f 


Not operational in current program version. 
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5. GUIDELINES FOR PROGRAM USE 


Developers of nonlinear programs have long cautioned prospective program 
users that the successful use of a large, nonlinear program is considerably 
different than the use of a linear program; in fact, it is often a tedious, 
frustrating and expense exercise often requiring repeated attempts at solving 
a particular problem. This is due in part, to the required nonlinear solution 
algorithms which are often slowly convergent, and somtimes divergent. Various 
nonlinear kinematic and material models often introduce instabilities which 
are sensitive to numerical roundoff, grid modeling, numerical integration 
orders, and the like. This, coupled with the unfamilarity of most program 
users with the inner workings and idiosyncrasies of a nonlinear program, leads 
to many nonlinear computer programs being shelved after a short time. This 
section presents some guidelines which hopefully will prove beneficial to 
less-experienced users of nonlinear programs. The following paragraphs are 
not intended to be a recipe for obtaining a good solution, rather they are a 
collection of assorted guidelines which should always be considered when 
setting up a problem or when difficulties arise. 

For geometrically nonlinear (large strain) problems, it often is the case 
that at small strain levels where the gross deformations appear to be relatively 
linear with load, the program may experience considerable difficulty in con- 
verging to a solution. This is generally due to the initial coupling between 
bending, shear and membrane strains which are extremely nonlinear at low strain 
levels even though may contribute insignificantly to the total gross deformation. 
This often requires relatively small load steps at the onset during which 
these nonlinearities develop, after which larger steps may be taken. 

Difficulties in obtaining convergence may also occur at the onset of 
plastic deformation if a large portion of the structure undergoes simultaneous 
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plastic deformation. A similar situation when the strucutre has significant 
plastic zones and it experiences sudden load reversals causing the structure 
to unload elastically. This latter case causes discontinuities in the 
response of the structure and consequent difficulties in obtaining convergence. 
In cycles load cases, it is mandatory that two or three relatively small load 
steps be taken after the load reversal, otherwise solution convergence will 
most likely not be achieved. 

The solution scheme in the program allows the user to specify how often 
(load step frequency) the nonlinear stiffness is updated and how often equi- 
librium iteration is performed. If the nonlinear incremental stiffness is 
updated (reformed) infrequently, this will provide a more economical solution 
but will tend to introduce mathematical instability. If equilibrium iteration 
is not performed at each step, this will tend to stabilize the solution but 
will make it less accurate due to accumulated error from step to step. If 
in doubt, equilibrium iteration should be performed at each step. However, 
on a step where load reversal begins it is often advisable not to perform 
equilibrium iteration. 

For 4 to 8 node, two-dimensional elements, a 2 x 2 Gaussian integration 
order for the stiffness is generally sufficient. For 16 node, three- 
dimensional elements, 2x2x2 itegration is generally adequate; however, 
for 20 node elements a 3 x 3 x 3 order is necessary (a lower order integration 
will often provide inconsistent results). 

Users should be discouraged from using the three-dimensional isoparametric 
element to model thin plate or shell behavior. In general, it has been 
demonstrated that the element gives a rather poor representation of thin shell 
behavior. In particular, as the thickness is reduced, the three-dimensional 
solution becomes excessively stiff when compared to a Kirchhoff-type thin 
shell or plate solution. 
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The substructure mesh generator option ,(MG0 = 1) which is included in 
the program should be used with caution. The mesh generator may generate 
incorrect nodal coordinates, element connectivity and/or boundary conditions 
under certain conditions. Consequently, the user is cautioned to check 
generated grids carefully when using the MGO = 1 option and, preferably 
is urged to use the MGO = 0 option of mesh input. 
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6. PROGRAM RESOURCE REQUIREMENTS, CAPACITY, AND LIMITATIONS 


The AGGIE program has been developed .on a Amdahl 470 V/6 computer using 
the IBM 370 Fortran IV compiler and associated software. The minimum memory 
requirement is approximately 400k bytes when overlaid as shown in Table 2 and 
compiled in double precision. The memory requirement can be reduced to about 
300k if certain modules are not required for a particular problem or if 
additional overlaying is done. The program requires 12 scratch disk files 
for the usual nonlinear problem (plus 1 other if displacement incrementation 
is used); see comnent cards in the MAIN program for description of the 
required files. 

As noted previously, the majority of all problem size - dependent arrays 
are incorporated into a single common block called AMAT. The size of this 
common block is defined' in the MAIN program. If problem size dictates it, 
this common block should be expanded for more efficient program execution. 

The variable MTOT (in MAIN) should be set equal to the size of common block. 

In addition, the size of stiffness matrix block size (LBLOCK) is defined in 
MAIN. This portion of the common block used for 1 block of the stiffness 
matrix should be set as large as practical for efficient execution and 
decreased I/O time (however, it should not exceed approximately 75% of MTOT). 

As explained earlier in the summary of the AGGIE program, the capacity 
of the program is, on a practical basis, "unlimited". If the complete model 
data will not fit into the array space allocated in common AMAT, the structure 
may be separated into element groups. In addition, the program uses an out- 
of-core equation solver (if necessary) to solve problems whose equilibrium 
equation will not fit entirely in the main core memory (with common AMAT). 
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7. PROGRAM RESTART 


A problem involving a large number of degrees-of- freedom, a large number 
of load steps, or exhibiting complex nonlinear behavior can often be analyzed 
most efficiently by using the restart capability -of the program. This generally 
requires running the program for several load steps and having the program write 
and save restart information which can then be used later to restart the program 
at that point and continue the analysis. This allows the user to analyze the 
output and determine the adequacy and accuracy of the solution at an intermediate 

point and then continue it if necessary without recomputing the prior solution 

\ 

history. 

To perform a restart run, a first job must have been run which executed a 
successfully for NSTE load^or time steps) and which wrote and saved the necessary 
restart files. The files which must be saved are 1, 4, 7, 8, 9, and 10. The 
restart job is then executed by supplying the same input as used for the original 
run except that on the Solution Control Card II. 1, the value of MODEX is set 
equal to 2, the Value of TSTART is set equal to the new starting value (which 
msut be the same as the load (time) value at which the restart information was 
written), and the value of NSTE is set equal to the number of additional steps 
desired. If deemed necessary, the Nonlinear Solution Variables on Card II. 2 
may be changed for the restart run. All loading information on Cards VII I— XI 
may be modified, except that it should be noted that only loading information 
subsequent to TSTART is pertinent. No changes may be made to the element mesh, 
boundary inditions, material models, type of element nonlinearity, etc. 
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8. EXAMPLE PROBLEMS 


The example problems which follow were chosen to demonstrate the various 
facets of the input required to run a particular problem, to acquaint the 
user with the type of problems that can be solved, and to provide documented 
results which may be used to check out the program upon its initial instal- 
lation at a different computer facility. In an attempt to establish some 
confidence level in the program, the results are compared to either closed- 
form solutions, other finite element solutions or experimental results. 

The particular problems presented here do not by any means represent 
the "typical" problem that can be solved nor the upper capacity or capability 
of the program. For example, the capacity of the program in terms of 
number of elements allowed is practically unlimited. Additional problems 
which have been solved with the program are presented in Ref. 8. 

All solutions presented here were obtained with the computer program 
compiled in IBM Fortran 6 in double precision and executed on the Amdahl 
470 V/6 located at the Data Processing Center of Texas A&M University. 

Using the overlay structure discussed previously, the program requires 
about 450K bytes of storage (excluding disk buffer space). All quoted com- 
puter execution times refer to CPU time only and do not include items such' 
as I/O time, spooling time, virtual memory paging time, etc. 

Example Problem I 

The first example problem involves the large, plastic deformation of 
a simply supported curcular plate composed of high strainhardening 2024-0 
aluminum. The plate under consideration is centrally loaded by a circular 
punch of radius 0.1875 inches to a maximum load of 4,000 pounds. The 
dimensions of the plate are given in Fig. 22 . The unaxial material prop- 
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Fig. 22 Deflection history, circular plate with punch. 


erties for this plate are given in Ref. 9 along with the plate test data 
and theoretical results obtained with the PLANS computer program. The 
material uniaxial stress-strain curve was approximated by five linear 
segments with nodal engineering stress, strain values as follows: 10,000 

psi , 0.000995; 20,000 psi , 0.0131; 28,000 psi , 0.04502; 36',000 psi , 0.118; 
and 48,000 psi, 0.36311. The value of Poisson's ratio was not specified 
in Ref. 9; hence, a vlaue of 0.325 was chosen for use in the present analysis. 

Using a 100 pound load increment (40 steps), the plate was modeled 
using eleven and twenty- two 8-node isoparametric continuum elements with 
the results of both idealizations being the same indicating a converged 
solution. The predicted deflections of the lower surface of the plate at 
radii of 0.0 inch and 1.01 inches are presented in Fig. 22 along with 
theoretical and experimental data from Ref. 9. 

The PLANS results shown in Fig. 22 were obtained by a finite element 
idealization of fourteen annular shell elements and using a 5.5 pound 
load increment. The fact that the annular elements of PLANS include bending 
effects but not the effects of transverse shear and normal stresses, where- 
as the elements of AGGIE are pure continuum element may account for the 
majority of the discrepency between the results of the two programs. 

The theoretical results of both PLANS and AGGIE were obtained using 
pure kinematic hardening; however, since the loading is monotonically 
increasing and each point in the structure will experience a near radial 
stress loading path, the results are the same as if either isotropic hardening 
or the mechanical sublayer models had been used. 

The necessary input data for the eleven element idealization with 
kinematic hardening is shown in Fig. 23. 

The CPU execution time for the eleven element idealization was ap- 
proximately 46 seconds. 
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Fig. 23 Input Data for Example Problem 1 



Example Problem 2 

The second problem involves the same simply supported circular plate 
considered in Example Problem 1, but subjected to a cyclic loading history 
of 1325 pounds, then -773 pounds, and finally 1325 pounds. Since the 
loading is not monotonic, the predicted structural response will depend 
on the particular hardening rule used in the analysis. Based on the results 
of the previous sample problem, the plate was modeled using nine 8-node 
isoparametric (axisymmetric) elements. In modeling the plate, it was 
assumed that the plug of material under the punch acts rigidly; hence, 
the material under the punch was removed and roller boundary conditions 
inserted at the punch boundary. The punch load was then approximated 
as a shear load distributed through the thickness at the punch boundary. 
Since the loading was not severe enough to introduce appreciable 
geometric nonlinearities, the analyses were run as materially nonlinear 
only. 

The structural response predicted by AGGIE using the mechanical sub- 
layer model is shown in Fig. 24 along with experimental data from Ref. 10. 
The material constants included a Poisson's ratio 0.325, Young's modulus 

r 

of 10.05 x 10 psi , and a piecewise linear uniaxial stress-strain curve 
with the following nodal stress, strain values: 8,000 psi, 0.000796; 

16250 psi, 0.00575; 25,000 psi, 0.02925; and 35,000 psi, 0.10. Compared 
to experiment, the results are qualitatively good. The required program 
input is shown in Fig. 25. 

Results for combined kinematic-isotropic hardening with a constant 
ratio y = 0.5 are shown in Fig. 26. The stress-strain curve was idealized 
slightly differently as: 8,000 psi, 0.000796; 12,000 psi, 0.0025; 

18,000 psi, 0.00925; and 35,000 psi, 0.10. The results in thts case are 
somewhat better than obtained with the mechanical sublayer model. The 
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MECHANICAL SUBLAYER MODEL 

Fig. 24 Center Deflection of Curcular Plate Subjected to Cyclic Punch Load - Mechanical Sublayer Model 
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Fig. 25 Input Data for Example Problem 2 - Mechanical Sublayer Model 
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Fig. 26 Center Deflection of Circular Plate Subjected to Cyclic Punch Load 
Kinematic - Isotropic Hardening Model with y ~ 0.5. 
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Fig. 27 Input Data for Example Problem 2 - Combined Isotropic-Kinematic Hardening Model with y-0.5 
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Fig. 28 Yield Surface Size vs. Equivalent Uniaxial Plastic Strain for 2024-0 Aluminum 
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Fig. 29 Center Deflect Hw of Circular Plate Subjected to Cyclic Punch Load - Combined 
Kinematric- Isotropic Hardening Model with Experimental Variable Ratio y 
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Fig. 30 Input Data for Example Problem 2 - Combined Kinematic-Isotropic Hardening Model 
with Experimental Variable Ratio Data 
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Fig. 31 Center Deflection of Circular Plate Subjected to Cyclic Punch Load - Isotropic Hardening 
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Fig. 32 Center Deflection of Circular Plate Subjected to Cyclic Punch Load - Kinematic Hardening 



required input data is shown in Fig. 27 . 

For the 2024-0 aluminum used in this test, a yield surface size vs. 
plastic strain curve has been obtained from experimental cyclic test data 
provided by Winter.^ This experimental yield surface curve along with 
several others are shown in Fig. 28. This data has been used in the variable 
ratio, combined kinematic-isotropic hardening model (TPLAS.EQ.3) to obtain 
the results shown in Fig. 29. The results are qualitatively good. Input 
data for this solution is shown in Fig. 30. 

For comparison purposes, solutions obtained with the pure isotropic 
and pure kinematic hardening models are shown in Figs. 31 and 32. The stress 
strain curve used was the same as used in the kinematic hardening model. The 
results are poor as would' be expected. Since isotropic hardening predicts 
a reverse ideal Bauschinger effect, the structural response after the 
initial loading of 1325 pounds is predominately elastic, resulting in 
poor agreement with the experimental data. For kinematic hardening, the 
slight strainhardening rate predicted after a reversal of loading results 
in a very soft tangent stiffness and a corresponding overestimate of 
the displacement increments. Additional results for this problem may be 
found in Ref. 12. 

The CPU execution time for each of these problems was approximately 
44 seconds. 

Example Problem 3 

An infinitely long, thick-walled cylindrical tube made of incompressible 
material, 4 inches internal radius and 8 inches external radius, and 
subjected to an internal pressure of at most 40 psi is the next example 
problem. The material is assumed to be incompressible with Mooney-Rivlin 
material constants of C^=25 psi and psi. The finite element model 
consists of four 8-node axisymmetric elements with roller boundary conditions 
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Fig. 33 Cylinder under Internal pressure, Axi symmetric analysis. 
Displacement of interior surface vs. internal pressure 
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Fig. 34 Cylinder under internal pressure, Axi symmetric analysis. 
Displacement and hydrostatic pressure profiles 
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Fig. 35 Cylinder under internal pressure, Axi symmetric analysis. 

Stress profiles (radial and axial) 



C71 

CO 



Fig. 36 Cylinder under internal pressure, Axisymmetric analysis. 
Stress profile (tangential) 
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Fig. 37 Input Data for Example Problem 3 



at the top and bottom surface to simulate the infinite length conditions 
(i.e., no displacement allowed in the direction of the cylinder axis). 

Fig. 33 shows the displacement response obtained for a linear solution 
and nonlinear (total Lagrangian) solutions with conservative and noncon- 
servative loading. The nonlinear solution with non conservative loading 
is in excellent agreement with three closed-form exact solution; however, 
the conservative loading solution is considerably in error. This would 
be expected since the inner radius of the cylinder expends from 4 inches 
to about 10 inches and, hence, the surface area over which the pressure is 
applied increases by a factor of almost 6. A two element (8-node) solution 
yields virtually the same solution. 

Displacement, hydrostatic pressure, and stress profiles through the 
thickness for a two element (4 and 8 node) model are shown in Figs. 34 
through 36. Excellent agreement with exact solutions are noted in all 
cases except for some of the two, 4 node element solutions. 

The input data used for the 4 element (8-node) idealization is shown 
in Fig. 37. The CPU execution time for this data set was approximately 
7 seconds. 

Additional details on the solution of various incompressible material 
problems may be found in Ref. 5. 

Example Problem 4 

This problem is identical to the previous problem except that a thin 
steel jacket has been added around the thick rubber cylinder. The steel 
jacket is 1/4 inch thick with material properties E=30 x 10 6 psi , v=0.3, 
and an assumed uniaxial stress-strain curve approximated by stress, strain 
points: 50,000 psi, .00167; 60,000 psi, 0.01; and'65,000 psi, 0.10. 
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Fig. 38 Radial Displacement History for Thick Rubber-Steel Cylinder 
with Internal Pressure 
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The steel jacket is modeled with two 8-node elements and the isotropic 
hardening plasticity model (TPLAS.EQ.O) . The modeling of the rubber 
remains the same as in Example Problem 3. 

Figure 38 shows the radial displacement of the inner surface, 
rubber-steel interface, and outer surface as a function of pressure. 

The diplacements for the interface and outer surface are indistinguishable 
on this scale plot. It should be noted that the. steel jacket undergoes 
plastic deformation at a pressureof 1,725 psi. Linear thin shell membrane 
theory applied to the thin steel jacket by itself predicts the onset of 
plastic deformation at approximately 1,540 psi. The deformations in this 
case are considerably smaller than in the previous example; hence, this 
solution was obtained using conservative pressure loading. 

Example Problem 5 

The next example considers the dynamic response of a spherical cap 
subjected to a uniform normal pressure loading applied as a step function. 
The geometry and material properties are shown in Fig. 40. The material 
is aluminum and is assumed to have a yield stress of 24,000 psi with a 
secondary modulus E^ of 0.21 x 10 psi. The finite element model consists 
of ten 8-node axisymmetric elements. The solution is assumed to be geo- 
metrically and materially nonlinear 'with plasticity accounted for with the 
mechanical sublayer model. 

Figure 41 shows a comparison of the apex deflection solutions obtained 
with AGGIE (incremental solution), DYNAPLAS (based on this shell elements) 
and a pseudo force solution.^ Considering the different approximations in 
each solution, the agreement is quite good. Figure 42 presents displacement 
vs. time solutions for time steps of 2.5 and 5 y seconds indicating a 
relatively converged solution with respect to time step. 
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ALL SOLUTIONS: 

T.L. WITH MECHANICAL SUBLAYER 
At = 5 /i- sec. 



Fig. 41 Comparison of different nonlinear program 
solutions for a spherical cap. 
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T.L. with Mechanical Sublayer 

Pseudo Force Solutions - Houbolt Integration, Linear Extrap. 
Incremental Solutions- Newmark Integration, Equilibrium Iter. 



Fig. 42 Effects of time step size on nonlinear 
solutions for a spherical cap. 
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Fig. 43 Input Data for Example Problem 5 



Example Problem 6 

The next problem to be presented is the elastic-perfectly plastic 
analysis of a steel pressure vessel, as shown in Fig. 44 along with the 
finite element idealization used by AGGIE I. The results for the 
vertical deflection of point "A" predicted by the program using 45 8-node 
axisymmetric elements are shown in Fig. 45 along with experimental and 
theoretical results reported in Ref. 15. The material constants used were 
a Young's modulus of 29.12 x 10 psi , a Poisson's ratio of 0.3, 

and a yield stress of 40,540 psi. The finite element results presented 

15 

by Nayak and Zienkiewicz were obtained using fifty- three axisymmetric 
isoparametric elements. Since the deflections are small, the analysis 
considered only material nonlinearity which was modeled by the isotropic 
hardening rule. Sixteen load steps were used to obtain the solution with 
a CPU execution time of 103 seconds. Similar results were obtained using 
only five load steps. The input for the sixteen step solution is shown 
in Fig. 46. 
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Fig. 44 Ax i symmetric Steel Pressure Vessel Subjected to Uniform 
Internal Pressure 
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Pressure, 

ksi 



Fig. 45 Vertical Deflection of Point A for Axi symmetric Steel Pressure Vessel 
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Fig. 46 Input Data for Example Problem 6 (Continued) 
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Fig. 46 Input Data for Example Problem 6 (Concluded) 
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